# Heavy vs. light flowers # Callin Switzer # 12 October 2016 # update 8 Dec 2017 # # Compares the bees' frequency when they are buzzing on # heavy (metal added) vs. light flowers # Note: all pores were glued shut
#install packages
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp", "plyr", "effects", "gamm4", "viridis")
ipak(packages)## Loading required package: ggplot2## Loading required package: reshape2## Loading required package: lme4## Loading required package: Matrix## Loading required package: sjPlot## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package## Learn more about sjPlot with 'browseVignettes("sjPlot")'.## Loading required package: multcomp## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser## Loading required package: plyr## Loading required package: effects## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Loading required package: gamm4## Loading required package: mgcv## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## This is mgcv 1.8-20. For overview type 'help("mgcv-package")'.## This is gamm4 0.2-5## Loading required package: viridis## Loading required package: viridisLite## ggplot2 reshape2 lme4 sjPlot multcomp plyr effects gamm4
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## viridis
## TRUE# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))
# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"
# print system info
print(paste("last run ", Sys.time()))## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'## [1] "last run 2017-12-14 18:37:52"print(R.version)## _
## platform x86_64-apple-darwin15.6.0
## arch x86_64
## os darwin15.6.0
## system x86_64, darwin15.6.0
## status
## major 3
## minor 4.2
## year 2017
## month 09
## day 28
## svn rev 73368
## language R
## version.string R version 3.4.2 (2017-09-28)
## nickname Short SummerRead in data and do some visualizations
new_df = read.csv(file.path(dataDir, "02_HeavyLight_cleaned.csv"))
new_df$hive = as.factor(new_df$hive)
new_df$treatment <- mapvalues(new_df$treatment, from = c("sham", "weighted"), to = c("Sham", "Weighted"))
head(new_df)## freq amp datetime rewNum
## 1 250 0.17945 2016_09_28__09_59_13_112 1
## 2 300 0.56657 2016_09_28__09_59_14_520 2
## 3 400 0.68371 2016_09_28__09_59_15_030 3
## 4 410 0.80118 2016_09_28__09_59_15_707 4
## 5 370 0.52666 2016_09_28__09_59_16_277 5
## 6 360 0.31541 2016_09_28__09_59_17_096 6
## accFile beeID hive IT
## 1 2016_09_28__09_59_13_112_220_450_test.txt 7 3 3.98
## 2 2016_09_28__09_59_14_520_220_450_test.txt 7 3 3.98
## 3 2016_09_28__09_59_15_030_220_450_test.txt 7 3 3.98
## 4 2016_09_28__09_59_15_707_220_450_test.txt 7 3 3.98
## 5 2016_09_28__09_59_16_277_220_450_test.txt 7 3 3.98
## 6 2016_09_28__09_59_17_096_220_450_test.txt 7 3 3.98
## accFileAndFolder
## 1 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 2 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 3 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 4 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 5 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 6 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## dateFrozenOrMarked treatment amp_acc
## 1 28-Sep-16 Sham 17.64503
## 2 28-Sep-16 Sham 55.70993
## 3 28-Sep-16 Sham 67.22812
## 4 28-Sep-16 Sham 78.77876
## 5 28-Sep-16 Sham 51.78564
## 6 28-Sep-16 Sham 31.01377# calculate average number of each treatment
ll = list(tapply(new_df$treatment, new_df$beeID, table))
dd = data.frame(unlist(ll))
dd$beeTrt = row.names(dd)
dd$beeID = sapply(dd$beeTrt, FUN = function(x) strsplit(x, "\\.")[[1]][1])
dd$trt = sapply(dd$beeTrt, FUN = function(x) strsplit(x, "\\.")[[1]][2])
# shows average number of buzzes per treatment
tapply(dd$unlist.ll., INDEX = dd$trt, mean)## Sham Weighted
## 29.69444 35.86111tapply(dd$unlist.ll., INDEX = dd$trt, sd) # sd for number of buzzes## Sham Weighted
## 9.603034 15.224328# show histograms of frequencies
hist(new_df$freq, breaks = seq(215, 455, by = 10))
# Calculate average for each bee and for each treatment
frqMeans <- as.data.frame(tapply(X = new_df$freq, INDEX = list(new_df$beeID, new_df$treatment), mean))
frqMeans$beeID <- row.names(frqMeans)
frqMeans## Sham Weighted beeID
## 1 344.3333 347.0833 1
## 10 331.2000 382.5000 10
## 11 411.7391 387.6923 11
## 12 344.8276 339.2500 12
## 14 324.0625 281.6000 14
## 15 367.5000 376.3158 15
## 16 357.0588 386.0606 16
## 17 314.0741 341.6000 17
## 18 366.4706 345.0000 18
## 19 375.0000 379.2308 19
## 2 298.5294 322.9167 2
## 20 359.8148 313.3333 20
## 21 331.3158 310.7143 21
## 22 399.4286 388.2857 22
## 23 348.2759 351.2821 23
## 24 287.7778 273.3333 24
## 25 345.0000 352.3684 25
## 26 369.6000 392.7273 26
## 27 366.6667 350.0000 27
## 28 364.4444 374.1935 28
## 29 291.2500 319.3333 29
## 3 399.0323 388.9655 3
## 30 368.0769 369.0323 30
## 31 273.5294 320.4167 31
## 32 325.6000 246.1538 32
## 4 356.6667 371.3636 4
## 5 370.7692 353.7143 5
## 6 284.2857 326.5714 6
## 7 366.9697 396.3415 7
## 8 333.8462 342.0833 8
## 9 344.4000 384.4444 9
## blue 319.1429 333.8462 blue
## gold 374.5098 365.2174 gold
## orange 309.6875 319.1000 orange
## pink 385.3333 367.5000 pink
## white 319.7674 326.9444 white# convert to long format for ggplot-ing
frqLong <- melt(frqMeans, id.vars = "beeID", measure.vars = c("Sham", "Weighted"),
variable.name = "trt", value.name = "frq")
nrow(frqLong)## [1] 72ggplot(frqLong, aes(x = trt, y = frq)) +
geom_boxplot() +
geom_line(aes(group = beeID))+
geom_point()
ggplot(frqLong, aes(x = trt, y = frq)) +
geom_boxplot() +
labs(y = "Buzz Frequency (Hz)", x = "Flower treatment")
#ggsave(file.path(figDir, "heavyLight_frq.pdf"), width = 5, height = 4)
# print some descriptive statistics
nrow(new_df)## [1] 2360unique(new_df$hive)## [1] 3 4
## Levels: 3 4# Calculate average amplitude for each bee and for each treatment
ampMeans <- as.data.frame(tapply(X = new_df$amp_acc, INDEX = list(new_df$beeID, new_df$treatment), mean))
ampMeans$beeID <- row.names(ampMeans)
ampMeans## Sham Weighted beeID
## 1 79.41488 49.082432 1
## 10 48.45538 30.837297 10
## 11 34.56381 22.888359 11
## 12 39.69695 36.219420 12
## 14 28.06241 22.151111 14
## 15 43.84567 39.297029 15
## 16 50.08044 36.654450 16
## 17 38.58280 34.215988 17
## 18 31.56744 23.471063 18
## 19 43.66984 24.653556 19
## 2 35.11669 34.075563 2
## 20 47.39991 35.394684 20
## 21 44.71363 34.498420 21
## 22 32.36932 24.702542 22
## 23 63.43800 46.047727 23
## 24 41.36935 22.975869 24
## 25 36.78398 20.052036 25
## 26 28.09219 22.162838 26
## 27 39.66976 30.902595 27
## 28 44.74322 39.945697 28
## 29 43.33800 22.381776 29
## 3 49.02683 56.246838 3
## 30 54.47406 29.037523 30
## 31 57.17575 30.517863 31
## 32 35.02128 5.009379 32
## 4 34.84784 37.962188 4
## 5 51.37301 32.967355 5
## 6 25.06068 36.736087 6
## 7 42.30789 47.528431 7
## 8 71.53457 28.270977 8
## 9 38.79524 29.125278 9
## blue 31.83149 19.595870 blue
## gold 34.72996 20.867171 gold
## orange 52.83226 46.007837 orange
## pink 30.62350 39.548159 pink
## white 22.55167 15.381023 whiteampLong <- melt(ampMeans, id.vars = "beeID", measure.vars = c("Sham", "Weighted"),
variable.name = "trt", value.name = "amp")
ggplot(ampLong, aes(x = trt, y = amp)) +
geom_boxplot() +
geom_point() +
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment")
ggplot(ampLong, aes(x = trt, y = amp)) +
geom_boxplot() +
geom_point() +
geom_line(aes(group = beeID))
ggplot(ampLong, aes(x = trt, y = amp)) +
geom_boxplot() +
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment")
#ggsave(file.path(figDir, 'heavyLightAmp.pdf'), width = 5, heigh = 4)
nrow(ampLong)## [1] 72Modeling frequency with GAMM
Use BIC to select model (decide what interactions and covariates to use)
Make sure that REML = FALSE when comparing BIC values
Start with the biggest model of interest, and then see what predictors can be removed
rr = new_df$treatment[1]
df1 = new_df[1, ]
for(ii in 1:(nrow(new_df) -1)){
if(rr == new_df$treatment[ii + 1]) next
else{
df1 = rbind(df1, new_df[ii+1, ])
rr = new_df$treatment[ii + 1]
}
}
df1 <- df1[df1$rewNum != 1, ]
df1$trtSwitch = df1$rewNum
colnames(df1)## [1] "freq" "amp" "datetime"
## [4] "rewNum" "accFile" "beeID"
## [7] "hive" "IT" "accFileAndFolder"
## [10] "dateFrozenOrMarked" "treatment" "amp_acc"
## [13] "trtSwitch"df2 <- df1[, c("beeID", "hive", "trtSwitch", "IT")]
df2## beeID hive trtSwitch IT
## 34 7 3 34 3.98
## 111 white 3 37 3.70
## 178 1 4 25 4.07
## 232 8 4 25 3.84
## 284 19 3 27 4.01
## 353 29 3 31 3.78
## 402 16 3 34 4.11
## 478 26 3 26 3.72
## 536 32 3 26 3.72
## 573 15 4 25 3.61
## 635 24 4 25 3.55
## 687 10 3 26 4.09
## 745 30 4 27 4.19
## 814 21 3 39 3.52
## 897 28 4 28 3.57
## 963 6 3 36 4.55
## 1002 orange 3 33 4.81
## 1129 17 4 28 4.25
## 1205 25 3 27 3.73
## 1267 31 3 25 3.76
## 1336 blue 3 36 3.91
## 1391 23 3 30 4.18
## 1459 12 3 30 4.10
## 1530 3 3 32 4.27
## 1594 22 4 36 3.65
## 1673 4 3 45 4.27
## 1710 9 3 26 4.17
## 1763 5 3 27 3.76
## 1832 2 4 35 4.42
## 1930 11 3 27 3.92
## 1986 20 3 34 4.28
## 2065 14 4 26 4.23
## 2127 pink 4 31 4.29
## 2204 27 3 34 3.51
## 2253 18 3 29 3.12
## 2310 gold 3 24 4.13sum(duplicated(df2$beeID))## [1] 0newdf2 <- merge(new_df, df2)
head(newdf2)## beeID hive IT freq amp datetime rewNum
## 1 1 4 4.07 310 0.65224 2016_09_27__10_52_32_696 1
## 2 1 4 4.07 330 0.47733 2016_09_27__10_52_41_293 2
## 3 1 4 4.07 240 0.55067 2016_09_27__10_52_42_247 3
## 4 1 4 4.07 290 0.28977 2016_09_27__10_52_42_974 4
## 5 1 4 4.07 360 0.47311 2016_09_27__10_52_51_809 5
## 6 1 4 4.07 360 0.45364 2016_09_27__10_53_00_124 6
## accFile
## 1 2016_09_27__10_52_32_696_220_450_test.txt
## 2 2016_09_27__10_52_41_293_220_450_test.txt
## 3 2016_09_27__10_52_42_247_220_450_test.txt
## 4 2016_09_27__10_52_42_974_220_450_test.txt
## 5 2016_09_27__10_52_51_809_220_450_test.txt
## 6 2016_09_27__10_53_00_124_220_450_test.txt
## accFileAndFolder
## 1 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 2 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 3 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 4 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 5 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 6 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## dateFrozenOrMarked treatment amp_acc trtSwitch
## 1 27-Sep-16 Weighted 64.13373 25
## 2 27-Sep-16 Weighted 46.93510 25
## 3 27-Sep-16 Weighted 54.14651 25
## 4 27-Sep-16 Weighted 28.49263 25
## 5 27-Sep-16 Weighted 46.52016 25
## 6 27-Sep-16 Weighted 44.60570 25newdf2$visitNum_centered <- newdf2$rewNum - newdf2$trtSwitch
# get treatment
trts <- (sapply(unique(newdf2$beeID), FUN = function(x){
startt = newdf2[newdf2$beeID == x & newdf2$rewNum == 1, "treatment" ]
return(as.character(startt))
}))
df4 <- data.frame(beeID = unique(newdf2$beeID), trts = trts)
newdf3 <- merge(newdf2, df4)
newdf3$trt2 <- as.factor(paste(newdf3$trts, "first"))
head(newdf3)## beeID hive IT freq amp datetime rewNum
## 1 1 4 4.07 310 0.65224 2016_09_27__10_52_32_696 1
## 2 1 4 4.07 330 0.47733 2016_09_27__10_52_41_293 2
## 3 1 4 4.07 240 0.55067 2016_09_27__10_52_42_247 3
## 4 1 4 4.07 290 0.28977 2016_09_27__10_52_42_974 4
## 5 1 4 4.07 360 0.47311 2016_09_27__10_52_51_809 5
## 6 1 4 4.07 360 0.45364 2016_09_27__10_53_00_124 6
## accFile
## 1 2016_09_27__10_52_32_696_220_450_test.txt
## 2 2016_09_27__10_52_41_293_220_450_test.txt
## 3 2016_09_27__10_52_42_247_220_450_test.txt
## 4 2016_09_27__10_52_42_974_220_450_test.txt
## 5 2016_09_27__10_52_51_809_220_450_test.txt
## 6 2016_09_27__10_53_00_124_220_450_test.txt
## accFileAndFolder
## 1 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 2 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 3 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 4 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 5 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 6 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## dateFrozenOrMarked treatment amp_acc trtSwitch visitNum_centered
## 1 27-Sep-16 Weighted 64.13373 25 -24
## 2 27-Sep-16 Weighted 46.93510 25 -23
## 3 27-Sep-16 Weighted 54.14651 25 -22
## 4 27-Sep-16 Weighted 28.49263 25 -21
## 5 27-Sep-16 Weighted 46.52016 25 -20
## 6 27-Sep-16 Weighted 44.60570 25 -19
## trts trt2
## 1 Weighted Weighted first
## 2 Weighted Weighted first
## 3 Weighted Weighted first
## 4 Weighted Weighted first
## 5 Weighted Weighted first
## 6 Weighted Weighted first# GAMM
# start with gamm so I can show change by visit number
g01 = gamm4(freq ~ s(visitNum_centered, by = trt2) + IT + hive + treatment + trt2, random = ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
summary(g01$gam) # Summary for paper ##
## Family: gaussian
## Link function: identity
##
## Formula:
## freq ~ s(visitNum_centered, by = trt2) + IT + hive + treatment +
## trt2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 373.350 64.000 5.834 6.17e-09 ***
## IT -5.075 15.811 -0.321 0.748
## hive4 -4.208 11.433 -0.368 0.713
## treatmentWeighted 1.919 3.431 0.559 0.576
## trt2Weighted first -10.607 10.771 -0.985 0.325
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum_centered):trt2Sham first 1 1 1.470 0.225
## s(visitNum_centered):trt2Weighted first 1 1 1.015 0.314
##
## R-sq.(adj) = -0.000401
## lmer.REML = 25097 Scale est. = 2354.9 n = 2360summary(g01$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 25096.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3193 -0.4698 0.1892 0.6934 2.3814
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 9.602e+02 3.099e+01
## Xr.0 s(visitNum_centered):trt2Weighted first 0.000e+00 0.000e+00
## Xr s(visitNum_centered):trt2Sham first 3.028e-11 5.503e-06
## Residual 2.355e+03 4.853e+01
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 373.350 64.000 5.834
## XIT -5.075 15.811 -0.321
## Xhive4 -4.208 11.433 -0.368
## XtreatmentWeighted 1.919 3.431 0.559
## Xtrt2Weighted first -10.607 10.771 -0.985
## Xs(visitNum_centered):trt2Sham firstFx1 2.285 1.885 1.212
## Xs(visitNum_centered):trt2Weighted firstFx1 2.569 2.550 1.007
##
## Correlation of Fixed Effects:
## X(Int) XIT Xhive4 XtrtmW Xtr2Wf X(N_):2Sf
## XIT -0.992
## Xhive4 -0.043 -0.011
## XtrtmntWght -0.041 0.011 0.003
## Xtrt2Wghtdf -0.237 0.165 -0.014 0.051
## X(N_):2SfF1 0.049 -0.028 -0.008 -0.707 -0.038
## X(N_):2WfF1 -0.032 0.012 -0.002 0.658 0.063 -0.466dev.off()## null device
## 1g02 = gamm4(freq ~ s(visitNum_centered, by = trt2) + hive + treatment + trt2, random = ~(1|beeID), data = newdf3)
anova(g01$mer, g02$mer) #g02 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g02$mer: NULL
## g01$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g02$mer 10 25152 25210 -12566 25132
## g01$mer 11 25154 25218 -12566 25132 0.1168 1 0.7325summary(g02$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 25104.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3173 -0.4693 0.1882 0.6954 2.3808
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 9.331e+02 3.055e+01
## Xr.0 s(visitNum_centered):trt2Weighted first 0.000e+00 0.000e+00
## Xr s(visitNum_centered):trt2Sham first 1.379e-10 1.174e-05
## Residual 2.355e+03 4.853e+01
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 352.967 7.981 44.23
## Xhive4 -4.248 11.276 -0.38
## XtreatmentWeighted 1.939 3.431 0.57
## Xtrt2Weighted first -10.030 10.478 -0.96
## Xs(visitNum_centered):trt2Sham firstFx1 2.262 1.884 1.20
## Xs(visitNum_centered):trt2Weighted firstFx1 2.583 2.550 1.01
##
## Correlation of Fixed Effects:
## X(Int) Xhive4 XtrtmW Xtr2Wf X(N_):2Sf
## Xhive4 -0.426
## XtrtmntWght -0.242 0.004
## Xtrt2Wghtdf -0.587 -0.012 0.051
## X(N_):2SfF1 0.171 -0.009 -0.707 -0.034
## X(N_):2WfF1 -0.158 -0.002 0.658 0.063 -0.466g03 = gamm4(freq ~ s(visitNum_centered, by = trt2) + hive + trt2, random = ~(1|beeID), data = newdf3)
anova(g02$mer, g03$mer) #g03 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g03$mer: NULL
## g02$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g03$mer 9 25151 25203 -12566 25133
## g02$mer 10 25152 25210 -12566 25132 0.328 1 0.5668summary(g03$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 25109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3127 -0.4673 0.1923 0.6924 2.3917
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 935.9 30.59
## Xr.0 s(visitNum_centered):trt2Weighted first 0.0 0.00
## Xr s(visitNum_centered):trt2Sham first 0.0 0.00
## Residual 2354.1 48.52
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 354.059 7.754 45.66
## Xhive4 -4.271 11.292 -0.38
## Xtrt2Weighted first -10.330 10.479 -0.99
## Xs(visitNum_centered):trt2Sham firstFx1 3.015 1.332 2.26
## Xs(visitNum_centered):trt2Weighted firstFx1 1.635 1.919 0.85
##
## Correlation of Fixed Effects:
## X(Int) Xhive4 Xtr2Wf X(N_):2Sf
## Xhive4 -0.438
## Xtrt2Wghtdf -0.593 -0.012
## X(N_):2SfF1 0.000 -0.008 0.003
## X(N_):2WfF1 0.002 -0.005 0.039 0.000g04 = gamm4(freq ~ s(visitNum_centered, by = trt2) +trt2, random = ~(1|beeID), data = newdf3)
anova(g03$mer, g04$mer) # g04 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g04$mer: NULL
## g03$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g04$mer 8 25149 25195 -12566 25133
## g03$mer 9 25151 25203 -12566 25133 0.1558 1 0.693summary(g04$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 25115.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3139 -0.4690 0.1912 0.6937 2.3917
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 9.112e+02 3.019e+01
## Xr.0 s(visitNum_centered):trt2Weighted first 1.513e-09 3.890e-05
## Xr s(visitNum_centered):trt2Sham first 3.322e-12 1.823e-06
## Residual 2.354e+03 4.852e+01
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 352.774 6.883 51.25
## Xtrt2Weighted first -10.376 10.345 -1.00
## Xs(visitNum_centered):trt2Sham firstFx1 3.008 1.332 2.26
## Xs(visitNum_centered):trt2Weighted firstFx1 1.631 1.919 0.85
##
## Correlation of Fixed Effects:
## X(Int) Xtr2Wf X(N_):2Sf
## Xtrt2Wghtdf -0.665
## X(N_):2SfF1 -0.004 0.003
## X(N_):2WfF1 0.000 0.039 0.000g05 = gamm4(freq ~ s(visitNum_centered) +trt2, random = ~(1|beeID), data = newdf3)
anova(g04$mer, g05$mer) # g05 slightly better (according to BIC)## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g05$mer: NULL
## g04$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g05$mer 6 25145 25180 -12567 25133
## g04$mer 8 25149 25195 -12566 25133 0.3437 2 0.8421summary(g05$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## freq ~ s(visitNum_centered) + trt2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 352.783 6.875 51.312 <2e-16 ***
## trt2Weighted first -10.188 10.329 -0.986 0.324
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum_centered) 1 1 5.478 0.0193 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.00388
## lmer.REML = 25120 Scale est. = 2353.5 n = 2360g06 = gamm4(freq ~ s(visitNum_centered), random = ~(1|beeID), data = newdf3)
anova(g05$mer, g06$mer) # g05 slightly better (according to BIC)## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g06$mer: NULL
## g05$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g06$mer 5 25144 25173 -12567 25134
## g05$mer 6 25145 25180 -12567 25133 1.0153 1 0.3136summary(g06$mer) # g06 Better## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 25127.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3036 -0.4740 0.1917 0.6901 2.3998
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 9.078e+02 3.013e+01
## Xr s(visitNum_centered) 9.414e-11 9.703e-06
## Residual 2.354e+03 4.851e+01
## Number of obs: 2360, groups: beeID, 36; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 348.268 5.127 67.93
## Xs(visitNum_centered)Fx1 2.587 1.093 2.37
##
## Correlation of Fixed Effects:
## X(Int)
## Xs(vstN_)F1 0.018g07 = gamm4(freq ~ 1, random = ~(1|beeID), data = newdf3)
anova(g06$mer, g07$mer) # g07 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g07$mer: NULL
## g06$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g07$mer 3 25146 25163 -12570 25140
## g06$mer 5 25144 25173 -12567 25134 5.584 2 0.0613 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# looks like GAMM is not necessary -- no predictors significantly affect frequencyModeling frequency with lmer
Do a similar thing to GAMM above (but make it more interpretable)
# fit lmer that is similar to gamm above
# interaction that may be important, based on domain knowledge
m0 = lmer(freq ~ I(scale(visitNum_centered)) * trt2 + IT + hive + I(scale(visitNum_centered^2)) + treatment + (1|beeID), data = newdf3, REML = FALSE)
summary(m0)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## freq ~ I(scale(visitNum_centered)) * trt2 + IT + hive + I(scale(visitNum_centered^2)) +
## treatment + (1 | beeID)
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 25152.3 25210.0 -12566.1 25132.3 2350
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3176 -0.4731 0.1881 0.6953 2.3769
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 848.9 29.14
## Residual 2351.9 48.50
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 373.6133 60.3485 6.191
## I(scale(visitNum_centered)) 2.0637 2.4466 0.844
## trt2Weighted first -10.5531 10.1598 -1.039
## IT -5.1688 14.9133 -0.347
## hive4 -4.1822 10.7784 -0.388
## I(scale(visitNum_centered^2)) 0.1987 1.5795 0.126
## treatmentWeighted 2.1164 3.6713 0.576
## I(scale(visitNum_centered)):trt2Weighted first 0.5702 4.2740 0.133
##
## Correlation of Fixed Effects:
## (Intr) I((N_) trt2Wf IT hive4 I((N_^ trtmnW
## I(scl(vN_)) 0.023
## trt2Wghtdfr -0.237 -0.049
## IT -0.992 0.002 0.164
## hive4 -0.042 -0.018 -0.013 -0.012
## I(sc(N_^2)) 0.026 -0.638 0.029 -0.038 0.018
## trtmntWghtd -0.031 -0.737 0.061 -0.003 0.010 0.357
## I((N_)):2Wf -0.031 -0.842 0.070 0.004 0.011 0.453 0.820# no interaction
m1 = update(m0, .~. - IT)
BIC(m0, m1) # stay with m1 -- no interaction## df BIC
## m0 10 25209.95
## m1 9 25202.31anova(m0, m1) # agrees with BIC## Data: newdf3
## Models:
## m1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + I(scale(visitNum_centered^2)) +
## m1: treatment + (1 | beeID) + I(scale(visitNum_centered)):trt2
## m0: freq ~ I(scale(visitNum_centered)) * trt2 + IT + hive + I(scale(visitNum_centered^2)) +
## m0: treatment + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 9 25150 25202 -12566 25132
## m0 10 25152 25210 -12566 25132 0.1199 1 0.7291summary(m1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## freq ~ I(scale(visitNum_centered)) + trt2 + hive + I(scale(visitNum_centered^2)) +
## treatment + (1 | beeID) + I(scale(visitNum_centered)):trt2
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 25150.4 25202.3 -12566.2 25132.4 2351
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3166 -0.4731 0.1875 0.6966 2.3778
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 852.1 29.19
## Residual 2351.8 48.50
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 352.8677 7.6941 45.86
## I(scale(visitNum_centered)) 2.0655 2.4466 0.84
## trt2Weighted first -9.9744 10.0391 -0.99
## hive4 -4.2258 10.7968 -0.39
## I(scale(visitNum_centered^2)) 0.1783 1.5784 0.11
## treatmentWeighted 2.1121 3.6713 0.58
## I(scale(visitNum_centered)):trt2Weighted first 0.5745 4.2740 0.13
##
## Correlation of Fixed Effects:
## (Intr) I((N_) trt2Wf hive4 I((N_^ trtmnW
## I(scl(vN_)) 0.196
## trt2Wghtdfr -0.586 -0.050
## hive4 -0.424 -0.018 -0.011
## I(sc(N_^2)) -0.093 -0.639 0.036 0.018
## trtmntWghtd -0.268 -0.737 0.062 0.010 0.358
## I((N_)):2Wf -0.218 -0.842 0.071 0.011 0.454 0.820m1.1 <- update(m1, .~. - I(scale(visitNum_centered^2)))
anova(m1.1, m1) # keep 1.1## Data: newdf3
## Models:
## m1.1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +
## m1.1: (1 | beeID) + I(scale(visitNum_centered)):trt2
## m1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + I(scale(visitNum_centered^2)) +
## m1: treatment + (1 | beeID) + I(scale(visitNum_centered)):trt2
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1.1 8 25148 25195 -12566 25132
## m1 9 25150 25202 -12566 25132 0.0128 1 0.9101summary(m1.1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +
## (1 | beeID) + I(scale(visitNum_centered)):trt2
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 25148.4 25194.6 -12566.2 25132.4 2352
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3159 -0.4706 0.1876 0.6974 2.3783
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 851.6 29.18
## Residual 2351.9 48.50
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 352.9488 7.6583 46.09
## I(scale(visitNum_centered)) 2.2418 1.8826 1.19
## trt2Weighted first -10.0149 10.0298 -1.00
## hive4 -4.2477 10.7920 -0.39
## treatmentWeighted 1.9640 3.4286 0.57
## I(scale(visitNum_centered)):trt2Weighted first 0.3557 3.8085 0.09
##
## Correlation of Fixed Effects:
## (Intr) I((N_) trt2Wf hive4 trtmnW
## I(scl(vN_)) 0.178
## trt2Wghtdfr -0.586 -0.035
## hive4 -0.425 -0.009 -0.012
## trtmntWghtd -0.252 -0.707 0.053 0.004
## I((N_)):2Wf -0.198 -0.806 0.061 0.003 0.790m2 <- update(m1.1, .~. - I(scale(visitNum_centered)):trt2)
anova(m1.1, m2) # keep m2 (based on BIC)## Data: newdf3
## Models:
## m2: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +
## m2: (1 | beeID)
## m1.1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +
## m1.1: (1 | beeID) + I(scale(visitNum_centered)):trt2
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 7 25146 25187 -12566 25132
## m1.1 8 25148 25195 -12566 25132 0.0087 1 0.9256summary(m2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +
## (1 | beeID)
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 25146.4 25186.8 -12566.2 25132.4 2353
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3167 -0.4685 0.1872 0.6957 2.3784
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 852.1 29.19
## Residual 2351.9 48.50
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 353.090 7.509 47.02
## I(scale(visitNum_centered)) 2.384 1.115 2.14
## trt2Weighted first -10.072 10.014 -1.01
## hive4 -4.251 10.795 -0.39
## treatmentWeighted 1.711 2.102 0.81
##
## Correlation of Fixed Effects:
## (Intr) I((N_) trt2Wf hive4
## I(scl(vN_)) 0.032
## trt2Wghtdfr -0.587 0.024
## hive4 -0.433 -0.011 -0.012
## trtmntWghtd -0.159 -0.194 0.007 0.002m3 <- update(m2, .~. - hive)
anova(m2, m3) # keep m3## Data: newdf3
## Models:
## m3: freq ~ I(scale(visitNum_centered)) + trt2 + treatment + (1 |
## m3: beeID)
## m2: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +
## m2: (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 25145 25179 -12566 25133
## m2 7 25146 25187 -12566 25132 0.1547 1 0.6941summary(m3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + treatment + (1 |
## beeID)
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 25144.6 25179.2 -12566.3 25132.6 2354
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3192 -0.4693 0.1879 0.6954 2.3798
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 855.9 29.26
## Residual 2351.9 48.50
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 351.811 6.784 51.86
## I(scale(visitNum_centered)) 2.379 1.115 2.13
## trt2Weighted first -10.121 10.035 -1.01
## treatmentWeighted 1.713 2.102 0.81
##
## Correlation of Fixed Effects:
## (Intr) I((N_) trt2Wf
## I(scl(vN_)) 0.031
## trt2Wghtdfr -0.657 0.024
## trtmntWghtd -0.175 -0.194 0.007m3.1 <- update(m3, .~. - treatment)
anova(m3, m3.1) # keep m3.1## Data: newdf3
## Models:
## m3.1: freq ~ I(scale(visitNum_centered)) + trt2 + (1 | beeID)
## m3: freq ~ I(scale(visitNum_centered)) + trt2 + treatment + (1 |
## m3: beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3.1 5 25143 25172 -12567 25133
## m3 6 25145 25179 -12566 25133 0.664 1 0.4151summary(m3.1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + (1 | beeID)
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 25143.2 25172.1 -12566.6 25133.2 2355
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2984 -0.4723 0.1895 0.6903 2.4024
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 856.2 29.26
## Residual 2352.5 48.50
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 352.780 6.680 52.81
## I(scale(visitNum_centered)) 2.556 1.094 2.34
## trt2Weighted first -10.181 10.036 -1.01
##
## Correlation of Fixed Effects:
## (Intr) I((N_)
## I(scl(vN_)) -0.003
## trt2Wghtdfr -0.666 0.025m4 <- update(m3.1, .~. - trt2)
summary(m4)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + (1 | beeID)
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 25142.3 25165.3 -12567.1 25134.3 2356
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3033 -0.4735 0.1920 0.6906 2.3989
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 881.3 29.69
## Residual 2352.5 48.50
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 348.269 5.054 68.91
## I(scale(visitNum_centered)) 2.586 1.093 2.36
##
## Correlation of Fixed Effects:
## (Intr)
## I(scl(vN_)) 0.018anova(m3.1, m4) #keep m4## Data: newdf3
## Models:
## m4: freq ~ I(scale(visitNum_centered)) + (1 | beeID)
## m3.1: freq ~ I(scale(visitNum_centered)) + trt2 + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 4 25142 25165 -12567 25134
## m3.1 5 25143 25172 -12567 25133 1.0153 1 0.3136m5 <- update(m4, .~. -I(scale(visitNum_centered)))
anova(m4, m5) # keep m5 -- baesd on BIC (agrees with GAMM analysis above)## Data: newdf3
## Models:
## m5: freq ~ (1 | beeID)
## m4: freq ~ I(scale(visitNum_centered)) + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 3 25146 25163 -12570 25140
## m4 4 25142 25165 -12567 25134 5.584 1 0.01813 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# final freq mod for paper
finmod <- update(m5, .~., REML = TRUE)
summary(finmod) # final model for paper -- same as gamm## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ (1 | beeID)
## Data: newdf3
##
## REML criterion at convergence: 25134.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3401 -0.4580 0.1806 0.6919 2.3616
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 901.6 30.03
## Residual 2358.4 48.56
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 348.056 5.109 68.12summary(finmod) # summary for paper -- no predictors predict frequency (when accounting for beeID)## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ (1 | beeID)
## Data: newdf3
##
## REML criterion at convergence: 25134.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3401 -0.4580 0.1806 0.6919 2.3616
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 901.6 30.03
## Residual 2358.4 48.56
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 348.056 5.109 68.12plot(finmod)
qqnorm(ranef(finmod)$beeID[[1]])
qqline(ranef(finmod)$beeID[[1]])
sjp.lmer(finmod, type = "re", sort = TRUE) # plot random effects to find any outliers## Plotting random effects...
# sjp.lmer(finmod, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effectsBootstrap CI’s for figure for paper
# set number of bootstrap replicates for models
nbootSims = 1000
table(new_df$hive) # more trials from hive 3##
## 3 4
## 1644 716# using hive 3, since it's the one with the most data
# however, hive doesn't affect model anyway
# calculate an average IT for prediction
ITmean = mean(tapply(new_df$IT, INDEX = new_df$beeID, FUN = function(x) x[1] ))
pframe <- expand.grid(trt2 = levels(droplevels(newdf3$trt2)),
IT = ITmean,
visitNum_centered = -50:50,
hive = factor(3, levels = levels(newdf3$hive)),
beeID = 99999,
treatment = levels(newdf3$treatment))
pframe$freq <- 0
pframe <- pframe[(pframe$trt2 == "Sham first" &
pframe$treatment == "Sham" &
pframe$visitNum_centered < 0) |
(pframe$trt2 == "Sham first" &
pframe$treatment == "Weighted" &
pframe$visitNum_centered > 0) |
(pframe$trt2 == "Weighted first" &
pframe$treatment == "Weighted" &
pframe$visitNum_centered < 0) |
(pframe$trt2 == "Weighted first" &
pframe$treatment == "Sham" &
pframe$visitNum_centered > 0)
, ]
# make a model to show the effect of treatment in a plot
m1 <- update(finmod, .~. + treatment)
pp <- predict(m1, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m1, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)## Warning in optwrap(object@optinfo$optimizer, ff, x0, lower = lower, control
## = control$optCtrl, : convergence code 3 from bobyqa: bobyqa -- a trust
## region step failed to reduce qprint(paste("Number of bootstrap samples", nrow(bb2$t)))## [1] "Number of bootstrap samples 1000"bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe## trt2 IT visitNum_centered hive beeID treatment freq
## 1 Sham first 3.965833 -50 3 99999 Sham 0
## 3 Sham first 3.965833 -49 3 99999 Sham 0
## 5 Sham first 3.965833 -48 3 99999 Sham 0
## 7 Sham first 3.965833 -47 3 99999 Sham 0
## 9 Sham first 3.965833 -46 3 99999 Sham 0
## 11 Sham first 3.965833 -45 3 99999 Sham 0
## 13 Sham first 3.965833 -44 3 99999 Sham 0
## 15 Sham first 3.965833 -43 3 99999 Sham 0
## 17 Sham first 3.965833 -42 3 99999 Sham 0
## 19 Sham first 3.965833 -41 3 99999 Sham 0
## 21 Sham first 3.965833 -40 3 99999 Sham 0
## 23 Sham first 3.965833 -39 3 99999 Sham 0
## 25 Sham first 3.965833 -38 3 99999 Sham 0
## 27 Sham first 3.965833 -37 3 99999 Sham 0
## 29 Sham first 3.965833 -36 3 99999 Sham 0
## 31 Sham first 3.965833 -35 3 99999 Sham 0
## 33 Sham first 3.965833 -34 3 99999 Sham 0
## 35 Sham first 3.965833 -33 3 99999 Sham 0
## 37 Sham first 3.965833 -32 3 99999 Sham 0
## 39 Sham first 3.965833 -31 3 99999 Sham 0
## 41 Sham first 3.965833 -30 3 99999 Sham 0
## 43 Sham first 3.965833 -29 3 99999 Sham 0
## 45 Sham first 3.965833 -28 3 99999 Sham 0
## 47 Sham first 3.965833 -27 3 99999 Sham 0
## 49 Sham first 3.965833 -26 3 99999 Sham 0
## 51 Sham first 3.965833 -25 3 99999 Sham 0
## 53 Sham first 3.965833 -24 3 99999 Sham 0
## 55 Sham first 3.965833 -23 3 99999 Sham 0
## 57 Sham first 3.965833 -22 3 99999 Sham 0
## 59 Sham first 3.965833 -21 3 99999 Sham 0
## 61 Sham first 3.965833 -20 3 99999 Sham 0
## 63 Sham first 3.965833 -19 3 99999 Sham 0
## 65 Sham first 3.965833 -18 3 99999 Sham 0
## 67 Sham first 3.965833 -17 3 99999 Sham 0
## 69 Sham first 3.965833 -16 3 99999 Sham 0
## 71 Sham first 3.965833 -15 3 99999 Sham 0
## 73 Sham first 3.965833 -14 3 99999 Sham 0
## 75 Sham first 3.965833 -13 3 99999 Sham 0
## 77 Sham first 3.965833 -12 3 99999 Sham 0
## 79 Sham first 3.965833 -11 3 99999 Sham 0
## 81 Sham first 3.965833 -10 3 99999 Sham 0
## 83 Sham first 3.965833 -9 3 99999 Sham 0
## 85 Sham first 3.965833 -8 3 99999 Sham 0
## 87 Sham first 3.965833 -7 3 99999 Sham 0
## 89 Sham first 3.965833 -6 3 99999 Sham 0
## 91 Sham first 3.965833 -5 3 99999 Sham 0
## 93 Sham first 3.965833 -4 3 99999 Sham 0
## 95 Sham first 3.965833 -3 3 99999 Sham 0
## 97 Sham first 3.965833 -2 3 99999 Sham 0
## 99 Sham first 3.965833 -1 3 99999 Sham 0
## 104 Weighted first 3.965833 1 3 99999 Sham 0
## 106 Weighted first 3.965833 2 3 99999 Sham 0
## 108 Weighted first 3.965833 3 3 99999 Sham 0
## 110 Weighted first 3.965833 4 3 99999 Sham 0
## 112 Weighted first 3.965833 5 3 99999 Sham 0
## 114 Weighted first 3.965833 6 3 99999 Sham 0
## 116 Weighted first 3.965833 7 3 99999 Sham 0
## 118 Weighted first 3.965833 8 3 99999 Sham 0
## 120 Weighted first 3.965833 9 3 99999 Sham 0
## 122 Weighted first 3.965833 10 3 99999 Sham 0
## 124 Weighted first 3.965833 11 3 99999 Sham 0
## 126 Weighted first 3.965833 12 3 99999 Sham 0
## 128 Weighted first 3.965833 13 3 99999 Sham 0
## 130 Weighted first 3.965833 14 3 99999 Sham 0
## 132 Weighted first 3.965833 15 3 99999 Sham 0
## 134 Weighted first 3.965833 16 3 99999 Sham 0
## 136 Weighted first 3.965833 17 3 99999 Sham 0
## 138 Weighted first 3.965833 18 3 99999 Sham 0
## 140 Weighted first 3.965833 19 3 99999 Sham 0
## 142 Weighted first 3.965833 20 3 99999 Sham 0
## 144 Weighted first 3.965833 21 3 99999 Sham 0
## 146 Weighted first 3.965833 22 3 99999 Sham 0
## 148 Weighted first 3.965833 23 3 99999 Sham 0
## 150 Weighted first 3.965833 24 3 99999 Sham 0
## 152 Weighted first 3.965833 25 3 99999 Sham 0
## 154 Weighted first 3.965833 26 3 99999 Sham 0
## 156 Weighted first 3.965833 27 3 99999 Sham 0
## 158 Weighted first 3.965833 28 3 99999 Sham 0
## 160 Weighted first 3.965833 29 3 99999 Sham 0
## 162 Weighted first 3.965833 30 3 99999 Sham 0
## 164 Weighted first 3.965833 31 3 99999 Sham 0
## 166 Weighted first 3.965833 32 3 99999 Sham 0
## 168 Weighted first 3.965833 33 3 99999 Sham 0
## 170 Weighted first 3.965833 34 3 99999 Sham 0
## 172 Weighted first 3.965833 35 3 99999 Sham 0
## 174 Weighted first 3.965833 36 3 99999 Sham 0
## 176 Weighted first 3.965833 37 3 99999 Sham 0
## 178 Weighted first 3.965833 38 3 99999 Sham 0
## 180 Weighted first 3.965833 39 3 99999 Sham 0
## 182 Weighted first 3.965833 40 3 99999 Sham 0
## 184 Weighted first 3.965833 41 3 99999 Sham 0
## 186 Weighted first 3.965833 42 3 99999 Sham 0
## 188 Weighted first 3.965833 43 3 99999 Sham 0
## 190 Weighted first 3.965833 44 3 99999 Sham 0
## 192 Weighted first 3.965833 45 3 99999 Sham 0
## 194 Weighted first 3.965833 46 3 99999 Sham 0
## 196 Weighted first 3.965833 47 3 99999 Sham 0
## 198 Weighted first 3.965833 48 3 99999 Sham 0
## 200 Weighted first 3.965833 49 3 99999 Sham 0
## 202 Weighted first 3.965833 50 3 99999 Sham 0
## 204 Weighted first 3.965833 -50 3 99999 Weighted 0
## 206 Weighted first 3.965833 -49 3 99999 Weighted 0
## 208 Weighted first 3.965833 -48 3 99999 Weighted 0
## 210 Weighted first 3.965833 -47 3 99999 Weighted 0
## 212 Weighted first 3.965833 -46 3 99999 Weighted 0
## 214 Weighted first 3.965833 -45 3 99999 Weighted 0
## 216 Weighted first 3.965833 -44 3 99999 Weighted 0
## 218 Weighted first 3.965833 -43 3 99999 Weighted 0
## 220 Weighted first 3.965833 -42 3 99999 Weighted 0
## 222 Weighted first 3.965833 -41 3 99999 Weighted 0
## 224 Weighted first 3.965833 -40 3 99999 Weighted 0
## 226 Weighted first 3.965833 -39 3 99999 Weighted 0
## 228 Weighted first 3.965833 -38 3 99999 Weighted 0
## 230 Weighted first 3.965833 -37 3 99999 Weighted 0
## 232 Weighted first 3.965833 -36 3 99999 Weighted 0
## 234 Weighted first 3.965833 -35 3 99999 Weighted 0
## 236 Weighted first 3.965833 -34 3 99999 Weighted 0
## 238 Weighted first 3.965833 -33 3 99999 Weighted 0
## 240 Weighted first 3.965833 -32 3 99999 Weighted 0
## 242 Weighted first 3.965833 -31 3 99999 Weighted 0
## 244 Weighted first 3.965833 -30 3 99999 Weighted 0
## 246 Weighted first 3.965833 -29 3 99999 Weighted 0
## 248 Weighted first 3.965833 -28 3 99999 Weighted 0
## 250 Weighted first 3.965833 -27 3 99999 Weighted 0
## 252 Weighted first 3.965833 -26 3 99999 Weighted 0
## 254 Weighted first 3.965833 -25 3 99999 Weighted 0
## 256 Weighted first 3.965833 -24 3 99999 Weighted 0
## 258 Weighted first 3.965833 -23 3 99999 Weighted 0
## 260 Weighted first 3.965833 -22 3 99999 Weighted 0
## 262 Weighted first 3.965833 -21 3 99999 Weighted 0
## 264 Weighted first 3.965833 -20 3 99999 Weighted 0
## 266 Weighted first 3.965833 -19 3 99999 Weighted 0
## 268 Weighted first 3.965833 -18 3 99999 Weighted 0
## 270 Weighted first 3.965833 -17 3 99999 Weighted 0
## 272 Weighted first 3.965833 -16 3 99999 Weighted 0
## 274 Weighted first 3.965833 -15 3 99999 Weighted 0
## 276 Weighted first 3.965833 -14 3 99999 Weighted 0
## 278 Weighted first 3.965833 -13 3 99999 Weighted 0
## 280 Weighted first 3.965833 -12 3 99999 Weighted 0
## 282 Weighted first 3.965833 -11 3 99999 Weighted 0
## 284 Weighted first 3.965833 -10 3 99999 Weighted 0
## 286 Weighted first 3.965833 -9 3 99999 Weighted 0
## 288 Weighted first 3.965833 -8 3 99999 Weighted 0
## 290 Weighted first 3.965833 -7 3 99999 Weighted 0
## 292 Weighted first 3.965833 -6 3 99999 Weighted 0
## 294 Weighted first 3.965833 -5 3 99999 Weighted 0
## 296 Weighted first 3.965833 -4 3 99999 Weighted 0
## 298 Weighted first 3.965833 -3 3 99999 Weighted 0
## 300 Weighted first 3.965833 -2 3 99999 Weighted 0
## 302 Weighted first 3.965833 -1 3 99999 Weighted 0
## 305 Sham first 3.965833 1 3 99999 Weighted 0
## 307 Sham first 3.965833 2 3 99999 Weighted 0
## 309 Sham first 3.965833 3 3 99999 Weighted 0
## 311 Sham first 3.965833 4 3 99999 Weighted 0
## 313 Sham first 3.965833 5 3 99999 Weighted 0
## 315 Sham first 3.965833 6 3 99999 Weighted 0
## 317 Sham first 3.965833 7 3 99999 Weighted 0
## 319 Sham first 3.965833 8 3 99999 Weighted 0
## 321 Sham first 3.965833 9 3 99999 Weighted 0
## 323 Sham first 3.965833 10 3 99999 Weighted 0
## 325 Sham first 3.965833 11 3 99999 Weighted 0
## 327 Sham first 3.965833 12 3 99999 Weighted 0
## 329 Sham first 3.965833 13 3 99999 Weighted 0
## 331 Sham first 3.965833 14 3 99999 Weighted 0
## 333 Sham first 3.965833 15 3 99999 Weighted 0
## 335 Sham first 3.965833 16 3 99999 Weighted 0
## 337 Sham first 3.965833 17 3 99999 Weighted 0
## 339 Sham first 3.965833 18 3 99999 Weighted 0
## 341 Sham first 3.965833 19 3 99999 Weighted 0
## 343 Sham first 3.965833 20 3 99999 Weighted 0
## 345 Sham first 3.965833 21 3 99999 Weighted 0
## 347 Sham first 3.965833 22 3 99999 Weighted 0
## 349 Sham first 3.965833 23 3 99999 Weighted 0
## 351 Sham first 3.965833 24 3 99999 Weighted 0
## 353 Sham first 3.965833 25 3 99999 Weighted 0
## 355 Sham first 3.965833 26 3 99999 Weighted 0
## 357 Sham first 3.965833 27 3 99999 Weighted 0
## 359 Sham first 3.965833 28 3 99999 Weighted 0
## 361 Sham first 3.965833 29 3 99999 Weighted 0
## 363 Sham first 3.965833 30 3 99999 Weighted 0
## 365 Sham first 3.965833 31 3 99999 Weighted 0
## 367 Sham first 3.965833 32 3 99999 Weighted 0
## 369 Sham first 3.965833 33 3 99999 Weighted 0
## 371 Sham first 3.965833 34 3 99999 Weighted 0
## 373 Sham first 3.965833 35 3 99999 Weighted 0
## 375 Sham first 3.965833 36 3 99999 Weighted 0
## 377 Sham first 3.965833 37 3 99999 Weighted 0
## 379 Sham first 3.965833 38 3 99999 Weighted 0
## 381 Sham first 3.965833 39 3 99999 Weighted 0
## 383 Sham first 3.965833 40 3 99999 Weighted 0
## 385 Sham first 3.965833 41 3 99999 Weighted 0
## 387 Sham first 3.965833 42 3 99999 Weighted 0
## 389 Sham first 3.965833 43 3 99999 Weighted 0
## 391 Sham first 3.965833 44 3 99999 Weighted 0
## 393 Sham first 3.965833 45 3 99999 Weighted 0
## 395 Sham first 3.965833 46 3 99999 Weighted 0
## 397 Sham first 3.965833 47 3 99999 Weighted 0
## 399 Sham first 3.965833 48 3 99999 Weighted 0
## 401 Sham first 3.965833 49 3 99999 Weighted 0
## 403 Sham first 3.965833 50 3 99999 Weighted 0
## blo bhi predMean
## 1 335.8801 357.9694 346.6421
## 3 335.8801 357.9694 346.6421
## 5 335.8801 357.9694 346.6421
## 7 335.8801 357.9694 346.6421
## 9 335.8801 357.9694 346.6421
## 11 335.8801 357.9694 346.6421
## 13 335.8801 357.9694 346.6421
## 15 335.8801 357.9694 346.6421
## 17 335.8801 357.9694 346.6421
## 19 335.8801 357.9694 346.6421
## 21 335.8801 357.9694 346.6421
## 23 335.8801 357.9694 346.6421
## 25 335.8801 357.9694 346.6421
## 27 335.8801 357.9694 346.6421
## 29 335.8801 357.9694 346.6421
## 31 335.8801 357.9694 346.6421
## 33 335.8801 357.9694 346.6421
## 35 335.8801 357.9694 346.6421
## 37 335.8801 357.9694 346.6421
## 39 335.8801 357.9694 346.6421
## 41 335.8801 357.9694 346.6421
## 43 335.8801 357.9694 346.6421
## 45 335.8801 357.9694 346.6421
## 47 335.8801 357.9694 346.6421
## 49 335.8801 357.9694 346.6421
## 51 335.8801 357.9694 346.6421
## 53 335.8801 357.9694 346.6421
## 55 335.8801 357.9694 346.6421
## 57 335.8801 357.9694 346.6421
## 59 335.8801 357.9694 346.6421
## 61 335.8801 357.9694 346.6421
## 63 335.8801 357.9694 346.6421
## 65 335.8801 357.9694 346.6421
## 67 335.8801 357.9694 346.6421
## 69 335.8801 357.9694 346.6421
## 71 335.8801 357.9694 346.6421
## 73 335.8801 357.9694 346.6421
## 75 335.8801 357.9694 346.6421
## 77 335.8801 357.9694 346.6421
## 79 335.8801 357.9694 346.6421
## 81 335.8801 357.9694 346.6421
## 83 335.8801 357.9694 346.6421
## 85 335.8801 357.9694 346.6421
## 87 335.8801 357.9694 346.6421
## 89 335.8801 357.9694 346.6421
## 91 335.8801 357.9694 346.6421
## 93 335.8801 357.9694 346.6421
## 95 335.8801 357.9694 346.6421
## 97 335.8801 357.9694 346.6421
## 99 335.8801 357.9694 346.6421
## 104 335.8801 357.9694 346.6421
## 106 335.8801 357.9694 346.6421
## 108 335.8801 357.9694 346.6421
## 110 335.8801 357.9694 346.6421
## 112 335.8801 357.9694 346.6421
## 114 335.8801 357.9694 346.6421
## 116 335.8801 357.9694 346.6421
## 118 335.8801 357.9694 346.6421
## 120 335.8801 357.9694 346.6421
## 122 335.8801 357.9694 346.6421
## 124 335.8801 357.9694 346.6421
## 126 335.8801 357.9694 346.6421
## 128 335.8801 357.9694 346.6421
## 130 335.8801 357.9694 346.6421
## 132 335.8801 357.9694 346.6421
## 134 335.8801 357.9694 346.6421
## 136 335.8801 357.9694 346.6421
## 138 335.8801 357.9694 346.6421
## 140 335.8801 357.9694 346.6421
## 142 335.8801 357.9694 346.6421
## 144 335.8801 357.9694 346.6421
## 146 335.8801 357.9694 346.6421
## 148 335.8801 357.9694 346.6421
## 150 335.8801 357.9694 346.6421
## 152 335.8801 357.9694 346.6421
## 154 335.8801 357.9694 346.6421
## 156 335.8801 357.9694 346.6421
## 158 335.8801 357.9694 346.6421
## 160 335.8801 357.9694 346.6421
## 162 335.8801 357.9694 346.6421
## 164 335.8801 357.9694 346.6421
## 166 335.8801 357.9694 346.6421
## 168 335.8801 357.9694 346.6421
## 170 335.8801 357.9694 346.6421
## 172 335.8801 357.9694 346.6421
## 174 335.8801 357.9694 346.6421
## 176 335.8801 357.9694 346.6421
## 178 335.8801 357.9694 346.6421
## 180 335.8801 357.9694 346.6421
## 182 335.8801 357.9694 346.6421
## 184 335.8801 357.9694 346.6421
## 186 335.8801 357.9694 346.6421
## 188 335.8801 357.9694 346.6421
## 190 335.8801 357.9694 346.6421
## 192 335.8801 357.9694 346.6421
## 194 335.8801 357.9694 346.6421
## 196 335.8801 357.9694 346.6421
## 198 335.8801 357.9694 346.6421
## 200 335.8801 357.9694 346.6421
## 202 335.8801 357.9694 346.6421
## 204 339.2423 359.9579 349.2520
## 206 339.2423 359.9579 349.2520
## 208 339.2423 359.9579 349.2520
## 210 339.2423 359.9579 349.2520
## 212 339.2423 359.9579 349.2520
## 214 339.2423 359.9579 349.2520
## 216 339.2423 359.9579 349.2520
## 218 339.2423 359.9579 349.2520
## 220 339.2423 359.9579 349.2520
## 222 339.2423 359.9579 349.2520
## 224 339.2423 359.9579 349.2520
## 226 339.2423 359.9579 349.2520
## 228 339.2423 359.9579 349.2520
## 230 339.2423 359.9579 349.2520
## 232 339.2423 359.9579 349.2520
## 234 339.2423 359.9579 349.2520
## 236 339.2423 359.9579 349.2520
## 238 339.2423 359.9579 349.2520
## 240 339.2423 359.9579 349.2520
## 242 339.2423 359.9579 349.2520
## 244 339.2423 359.9579 349.2520
## 246 339.2423 359.9579 349.2520
## 248 339.2423 359.9579 349.2520
## 250 339.2423 359.9579 349.2520
## 252 339.2423 359.9579 349.2520
## 254 339.2423 359.9579 349.2520
## 256 339.2423 359.9579 349.2520
## 258 339.2423 359.9579 349.2520
## 260 339.2423 359.9579 349.2520
## 262 339.2423 359.9579 349.2520
## 264 339.2423 359.9579 349.2520
## 266 339.2423 359.9579 349.2520
## 268 339.2423 359.9579 349.2520
## 270 339.2423 359.9579 349.2520
## 272 339.2423 359.9579 349.2520
## 274 339.2423 359.9579 349.2520
## 276 339.2423 359.9579 349.2520
## 278 339.2423 359.9579 349.2520
## 280 339.2423 359.9579 349.2520
## 282 339.2423 359.9579 349.2520
## 284 339.2423 359.9579 349.2520
## 286 339.2423 359.9579 349.2520
## 288 339.2423 359.9579 349.2520
## 290 339.2423 359.9579 349.2520
## 292 339.2423 359.9579 349.2520
## 294 339.2423 359.9579 349.2520
## 296 339.2423 359.9579 349.2520
## 298 339.2423 359.9579 349.2520
## 300 339.2423 359.9579 349.2520
## 302 339.2423 359.9579 349.2520
## 305 339.2423 359.9579 349.2520
## 307 339.2423 359.9579 349.2520
## 309 339.2423 359.9579 349.2520
## 311 339.2423 359.9579 349.2520
## 313 339.2423 359.9579 349.2520
## 315 339.2423 359.9579 349.2520
## 317 339.2423 359.9579 349.2520
## 319 339.2423 359.9579 349.2520
## 321 339.2423 359.9579 349.2520
## 323 339.2423 359.9579 349.2520
## 325 339.2423 359.9579 349.2520
## 327 339.2423 359.9579 349.2520
## 329 339.2423 359.9579 349.2520
## 331 339.2423 359.9579 349.2520
## 333 339.2423 359.9579 349.2520
## 335 339.2423 359.9579 349.2520
## 337 339.2423 359.9579 349.2520
## 339 339.2423 359.9579 349.2520
## 341 339.2423 359.9579 349.2520
## 343 339.2423 359.9579 349.2520
## 345 339.2423 359.9579 349.2520
## 347 339.2423 359.9579 349.2520
## 349 339.2423 359.9579 349.2520
## 351 339.2423 359.9579 349.2520
## 353 339.2423 359.9579 349.2520
## 355 339.2423 359.9579 349.2520
## 357 339.2423 359.9579 349.2520
## 359 339.2423 359.9579 349.2520
## 361 339.2423 359.9579 349.2520
## 363 339.2423 359.9579 349.2520
## 365 339.2423 359.9579 349.2520
## 367 339.2423 359.9579 349.2520
## 369 339.2423 359.9579 349.2520
## 371 339.2423 359.9579 349.2520
## 373 339.2423 359.9579 349.2520
## 375 339.2423 359.9579 349.2520
## 377 339.2423 359.9579 349.2520
## 379 339.2423 359.9579 349.2520
## 381 339.2423 359.9579 349.2520
## 383 339.2423 359.9579 349.2520
## 385 339.2423 359.9579 349.2520
## 387 339.2423 359.9579 349.2520
## 389 339.2423 359.9579 349.2520
## 391 339.2423 359.9579 349.2520
## 393 339.2423 359.9579 349.2520
## 395 339.2423 359.9579 349.2520
## 397 339.2423 359.9579 349.2520
## 399 339.2423 359.9579 349.2520
## 401 339.2423 359.9579 349.2520
## 403 339.2423 359.9579 349.2520### Make frequency plots for paper
# predictions are from this model
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ (1 | beeID) + treatment
## Data: newdf3
##
## REML criterion at convergence: 25129.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3673 -0.4740 0.1874 0.6837 2.3307
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 901.5 30.02
## Residual 2357.8 48.56
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 346.642 5.230 66.28
## treatmentWeighted 2.610 2.064 1.26
##
## Correlation of Fixed Effects:
## (Intr)
## trtmntWghtd -0.214anova(m1, update(m1, .~. - treatment)) #p-val and BIC for paper## refitting model(s) with ML (instead of REML)## Data: newdf3
## Models:
## update(m1, . ~ . - treatment): freq ~ (1 | beeID)
## m1: freq ~ (1 | beeID) + treatment
## Df AIC BIC logLik deviance Chisq Chi Df
## update(m1, . ~ . - treatment) 3 25146 25163 -12570 25140
## m1 4 25146 25169 -12569 25138 1.5992 1
## Pr(>Chisq)
## update(m1, . ~ . - treatment)
## m1 0.206pframe$trt3 <- mapvalues(pframe$trt2, from = c("Sham first", "Weighted first"),
to = c("Sham -> Weighted", "Weighted -> Sham"))
g0 <- ggplot(pframe, aes(x=visitNum_centered, y=predMean))+
geom_line(aes(color = treatment))+
labs(y = "Frequency (Hz)", x = "Sonication number\n (0 is when treatment switched)") +
geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi, fill= treatment), alpha = 0.2) +
facet_wrap(~trt3) +
scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
scale_fill_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
theme(legend.position = c(0.5, 0.90),
legend.background = element_rect(fill=alpha('gray95', 1)),
legend.direction="horizontal") +
ylim(c(330, 370))
g0
ggsave(filename = file.path(figDir, "WeightedSham_Freq_Timeseries.pdf"), width =5, height = 3.6)
g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
aes(x = visitNum_centered, y = freq), position = position_jitter(height = 5), size = 0.1, pch = 16, alpha = 0.4) +
ylim(c(220, 470))## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.g1 # note that the distribution of frequency is a bit bimodal -- this could be 
# a problem, so we ran the analysis with a freq cutoff at 270
# we found no significant differences in the models either way.
ggsave(filename = file.path(figDir, "WeightedSham_Freq_Timeseries_RawData.pdf"), width =5, height = 3.6)
g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
aes(x = visitNum_centered, y = freq), position = position_jitter(height = 5), size = 0.1, pch = 16, alpha = 0.4) +
ylim(c(270, 480))## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.g1## Warning: Removed 730 rows containing missing values (geom_point).
ggsave(filename = file.path(figDir, "WeightedSham_Freq_Timeseries_RawData_trimmed.pdf"), width =5, height = 3.6)## Warning: Removed 738 rows containing missing values (geom_point).# mean and 95% CI (without visit number)
g0 <- ggplot(pframe[pframe$visitNum_centered %in% c(-25, 25), ], aes(x=treatment, y=predMean))+
geom_point(aes(color = treatment), alpha = 1) +
labs(y = "Frequency (Hz)", x = "Flower treatment") +
geom_errorbar(aes(x = treatment, ymin = blo, ymax = bhi, color= treatment), alpha = 1, width = 0.1) +
scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
theme(legend.position = "none",
legend.background = element_rect(fill=alpha('gray95', 1)),
legend.direction="horizontal")
g0
ggsave(filename = file.path(figDir, "WeightedSham_Freq_pointWhiskers.pdf"), width =5, height = 3.6)
## add raw data
g1 <- g0 + geom_point(data = newdf3, aes(x = treatment, y = freq, color = treatment), position = position_jitter(width = 0.1, height = 6), alpha = 0.4, pch = 16, size = 0.2)
g1
ggsave(filename = file.path(figDir, "WeightedSham_Freq_pointWhiskers_rawData.pdf"), width =5, height = 3.6)
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
g0 <- ggplot(pframe, aes(x=visitNum_centered, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "buzzNumber") +
geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi), alpha = 0.2) +
facet_wrap(~trt2)
g0
# plot effect to see if it agrees
plot(Effect(c("treatment"), m1))
Modeling amplitude with gamm
g01 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2) + trt2 + IT + hive + treatment, random = ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
summary(g01$gam) # Summary for paper ##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(amp_acc) ~ s(visitNum_centered, by = trt2) + trt2 + IT +
## hive + treatment
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.72977 0.58219 4.689 2.90e-06 ***
## trt2Weighted first -0.05297 0.09833 -0.539 0.590
## IT 0.21347 0.14365 1.486 0.137
## hive4 0.09955 0.10391 0.958 0.338
## treatmentWeighted -0.32857 0.04921 -6.678 3.02e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum_centered):trt2Sham first 1.414 1.414 4.924 0.0175 *
## s(visitNum_centered):trt2Weighted first 2.800 2.800 3.683 0.0323 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0571
## lmer.REML = 4557.5 Scale est. = 0.384 n = 2360summary(g01$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 4557.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9897 -0.5342 0.1125 0.6707 2.6828
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.076179 0.27601
## Xr.0 s(visitNum_centered):trt2Weighted first 0.376254 0.61340
## Xr s(visitNum_centered):trt2Sham first 0.008701 0.09328
## Residual 0.384002 0.61968
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 2.729772 0.582186 4.689
## Xtrt2Weighted first -0.052966 0.098331 -0.539
## XIT 0.213470 0.143651 1.486
## Xhive4 0.099547 0.103913 0.958
## XtreatmentWeighted -0.328574 0.049206 -6.678
## Xs(visitNum_centered):trt2Sham firstFx1 0.059922 0.036356 1.648
## Xs(visitNum_centered):trt2Weighted firstFx1 0.003639 0.188745 0.019
##
## Correlation of Fixed Effects:
## X(Int) Xtr2Wf XIT Xhive4 XtrtmW X(N_):2Sf
## Xtrt2Wghtdf -0.238
## XIT -0.991 0.166
## Xhive4 -0.041 -0.011 -0.012
## XtrtmntWght -0.062 0.061 0.015 0.000
## X(N_):2SfF1 0.051 -0.029 -0.030 -0.003 -0.455
## X(N_):2WfF1 -0.001 0.060 0.000 0.001 0.023 -0.011dev.off()## null device
## 1g011 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2) + IT + hive + treatment, random = ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g011$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
summary(g011$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(amp_acc) ~ s(visitNum_centered, by = trt2) + IT + hive +
## treatment
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.65539 0.55974 4.744 2.22e-06 ***
## IT 0.22614 0.14023 1.613 0.107
## hive4 0.09890 0.10286 0.961 0.336
## treatmentWeighted -0.32629 0.04905 -6.652 3.58e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum_centered):trt2Sham first 1.396 1.396 4.868 0.0185 *
## s(visitNum_centered):trt2Weighted first 2.824 2.824 3.711 0.0286 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0543
## lmer.REML = 4555 Scale est. = 0.38399 n = 2360summary(g011$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 4555
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9904 -0.5310 0.1120 0.6721 2.6775
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.074527 0.27300
## Xr.0 s(visitNum_centered):trt2Weighted first 0.387226 0.62227
## Xr s(visitNum_centered):trt2Sham first 0.008222 0.09068
## Residual 0.383989 0.61967
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 2.65539 0.55974 4.744
## XIT 0.22614 0.14023 1.613
## Xhive4 0.09890 0.10286 0.961
## XtreatmentWeighted -0.32629 0.04905 -6.652
## Xs(visitNum_centered):trt2Sham firstFx1 0.05943 0.03583 1.659
## Xs(visitNum_centered):trt2Weighted firstFx1 0.01018 0.19078 0.053
##
## Correlation of Fixed Effects:
## X(Int) XIT Xhive4 XtrtmW X(N_):2Sf
## XIT -0.994
## Xhive4 -0.045 -0.011
## XtrtmntWght -0.049 0.005 0.001
## X(N_):2SfF1 0.046 -0.026 -0.003 -0.463
## X(N_):2WfF1 0.014 -0.010 0.002 0.018 -0.008dev.off()## null device
## 1anova(g01$mer, g011$mer) # go with g011## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g011$mer: NULL
## g01$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g011$mer 10 4553.4 4611.1 -2266.7 4533.4
## g01$mer 11 4555.1 4618.5 -2266.5 4533.1 0.3497 1 0.5543g02 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2) + IT + treatment, random = ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g02$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
summary(g02$gam) # Summary for paper ##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(amp_acc) ~ s(visitNum_centered, by = trt2) + IT + treatment
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.67960 0.55806 4.802 1.67e-06 ***
## IT 0.22770 0.13995 1.627 0.104
## treatmentWeighted -0.32686 0.04912 -6.655 3.52e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum_centered):trt2Sham first 1.412 1.412 4.868 0.0182 *
## s(visitNum_centered):trt2Weighted first 2.810 2.810 3.669 0.0310 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0534
## lmer.REML = 4553.2 Scale est. = 0.384 n = 2360summary(g02$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 4553.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9949 -0.5328 0.1143 0.6739 2.6875
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.074211 0.27242
## Xr.0 s(visitNum_centered):trt2Weighted first 0.379482 0.61602
## Xr s(visitNum_centered):trt2Sham first 0.008663 0.09308
## Residual 0.384001 0.61968
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 2.679603 0.558065 4.802
## XIT 0.227704 0.139951 1.627
## XtreatmentWeighted -0.326857 0.049116 -6.655
## Xs(visitNum_centered):trt2Sham firstFx1 0.059463 0.036299 1.638
## Xs(visitNum_centered):trt2Weighted firstFx1 0.009651 0.189097 0.051
##
## Correlation of Fixed Effects:
## X(Int) XIT XtrtmW X(N_):2Sf
## XIT -0.995
## XtrtmntWght -0.049 0.005
## X(N_):2SfF1 0.046 -0.026 -0.455
## X(N_):2WfF1 0.014 -0.010 0.019 -0.009dev.off()## null device
## 1anova(g011$mer, g02$mer) # g02 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g02$mer: NULL
## g011$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g02$mer 9 4552.4 4604.3 -2267.2 4534.4
## g011$mer 10 4553.4 4611.1 -2266.7 4533.4 0.9778 1 0.3228g03 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2) + treatment, random = ~(1|beeID), data = newdf3)
anova(g02$mer, g03$mer) #keep g03## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g03$mer: NULL
## g02$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g03$mer 8 4553.1 4599.2 -2268.5 4537.1
## g02$mer 9 4552.4 4604.3 -2267.2 4534.4 2.691 1 0.1009summary(g03$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 4553.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9839 -0.5287 0.1117 0.6716 2.6967
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.077961 0.27921
## Xr.0 s(visitNum_centered):trt2Weighted first 0.397553 0.63052
## Xr s(visitNum_centered):trt2Sham first 0.005922 0.07696
## Residual 0.384011 0.61969
## Number of obs: 2360, groups: beeID, 36; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 3.58211 0.05438 65.87
## XtreatmentWeighted -0.32459 0.04854 -6.69
## Xs(visitNum_centered):trt2Sham firstFx1 0.06166 0.03320 1.86
## Xs(visitNum_centered):trt2Weighted firstFx1 0.01321 0.19301 0.07
##
## Correlation of Fixed Effects:
## X(Int) XtrtmW X(N_):2Sf
## XtrtmntWght -0.451
## X(N_):2SfF1 0.228 -0.510
## X(N_):2WfF1 0.040 0.017 -0.009g04 <- gamm4(log(amp_acc) ~ s(visitNum_centered) + treatment, random = ~(1|beeID), data = newdf3)
anova(g03$mer, g04$mer) #g04 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g04$mer: NULL
## g03$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g04$mer 6 4551.2 4585.8 -2269.6 4539.2
## g03$mer 8 4553.1 4599.2 -2268.5 4537.1 2.1425 2 0.3426summary(g04$gam) # approx p-value for smooth##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(amp_acc) ~ s(visitNum_centered) + treatment
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.56539 0.05047 70.64 <2e-16 ***
## treatmentWeighted -0.29794 0.02688 -11.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum_centered) 2.056 2.056 6.862 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0378
## lmer.REML = 4555.5 Scale est. = 0.38524 n = 2360summary(g04$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 4555.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0433 -0.5427 0.1059 0.6832 2.7605
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.07754 0.2785
## Xr s(visitNum_centered) 0.02027 0.1424
## Residual 0.38524 0.6207
## Number of obs: 2360, groups: beeID, 36; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 3.56539 0.05047 70.64
## XtreatmentWeighted -0.29794 0.02688 -11.08
## Xs(visitNum_centered)Fx1 0.02517 0.04290 0.59
##
## Correlation of Fixed Effects:
## X(Int) XtrtmW
## XtrtmntWght -0.293
## Xs(vstN_)F1 0.024 -0.053g44 <- gamm4(log(amp_acc) ~ visitNum_centered + treatment, random = ~(1|beeID), data = newdf3)
g05 <- gamm4(log(amp_acc) ~ treatment, random = ~(1|beeID), data = newdf3)
anova(g44$mer, g05$mer) #g44 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g05$mer: NULL
## g44$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g05$mer 4 4560.8 4583.9 -2276.4 4552.8
## g44$mer 5 4549.2 4578.0 -2269.6 4539.2 13.63 1 0.0002226 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1g06 <- gamm4(log(amp_acc) ~ visitNum_centered , random = ~(1|beeID), data = newdf3)
anova(g44$mer, g06$mer) #g044 better## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## g06$mer: NULL
## g44$mer: NULL
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g06$mer 4 4666.2 4689.3 -2329.1 4658.2
## g44$mer 5 4549.2 4578.0 -2269.6 4539.2 119.02 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# final model -- only treatment matters
ampMod <- g04
summary(ampMod$mer)## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 4555.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0433 -0.5427 0.1059 0.6832 2.7605
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.07754 0.2785
## Xr s(visitNum_centered) 0.02027 0.1424
## Residual 0.38524 0.6207
## Number of obs: 2360, groups: beeID, 36; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 3.56539 0.05047 70.64
## XtreatmentWeighted -0.29794 0.02688 -11.08
## Xs(visitNum_centered)Fx1 0.02517 0.04290 0.59
##
## Correlation of Fixed Effects:
## X(Int) XtrtmW
## XtrtmntWght -0.293
## Xs(vstN_)F1 0.024 -0.053summary(ampMod$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(amp_acc) ~ s(visitNum_centered) + treatment
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.56539 0.05047 70.64 <2e-16 ***
## treatmentWeighted -0.29794 0.02688 -11.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(visitNum_centered) 2.056 2.056 6.862 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0378
## lmer.REML = 4555.5 Scale est. = 0.38524 n = 2360plot(ampMod$mer)
par(mfrow = c(2,1))
plot(ampMod$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
dev.off()## null device
## 1Modeling amplitude with lmer
# log transform for acceleration helps model fit better
# start with large model (including interaction)
ma0 = lmer(log(amp_acc) ~ visitNum_centered*trt2 + hive +IT * treatment+ (1|beeID), data = newdf3, REML = FALSE)
summary(ma0)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## log(amp_acc) ~ visitNum_centered * trt2 + hive + IT * treatment +
## (1 | beeID)
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 4541.7 4599.3 -2260.8 4521.7 2350
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0318 -0.5397 0.1116 0.6862 2.6657
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.06647 0.2578
## Residual 0.38294 0.6188
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3181167 0.5736148 5.785
## visitNum_centered 0.0017926 0.0010587 1.693
## trt2Weighted first -0.0608425 0.0918382 -0.662
## hive4 0.0948197 0.0975308 0.972
## IT 0.0589548 0.1420758 0.415
## treatmentWeighted -1.4174293 0.3083754 -4.596
## visitNum_centered:trt2Weighted first -0.0004062 0.0020914 -0.194
## IT:treatmentWeighted 0.2808255 0.0780233 3.599
##
## Correlation of Fixed Effects:
## (Intr) vstNm_ trt2Wf hive4 IT trtmnW vN_:2f
## vstNm_cntrd -0.015
## trt2Wghtdfr -0.229 0.040
## hive4 -0.038 -0.014 -0.013
## IT -0.992 0.038 0.162 -0.014
## trtmntWghtd -0.312 0.142 0.014 -0.009 0.316
## vstNm_c:2Wf -0.022 -0.805 -0.027 0.005 -0.005 -0.002
## IT:trtmntWg 0.307 -0.241 -0.016 0.010 -0.317 -0.990 0.115ma1 = update(ma0, .~. - visitNum_centered:trt2)
BIC(ma0, ma1) # interaction, ma1, is better## df BIC
## ma0 10 4599.347
## ma1 9 4591.618# pval for interaction (in case paper needs it)
anova(ma0, ma1)## Data: newdf3
## Models:
## ma1: log(amp_acc) ~ visitNum_centered + trt2 + hive + IT + treatment +
## ma1: (1 | beeID) + IT:treatment
## ma0: log(amp_acc) ~ visitNum_centered * trt2 + hive + IT * treatment +
## ma0: (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma1 9 4539.7 4591.6 -2260.9 4521.7
## ma0 10 4541.7 4599.3 -2260.8 4521.7 0.0377 1 0.846summary(ma1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## log(amp_acc) ~ visitNum_centered + trt2 + hive + IT + treatment +
## (1 | beeID) + IT:treatment
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 4539.7 4591.6 -2260.9 4521.7 2351
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0255 -0.5401 0.1127 0.6850 2.6689
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.06652 0.2579
## Residual 0.38294 0.6188
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3156701 0.5736617 5.780
## visitNum_centered 0.0016270 0.0006282 2.590
## trt2Weighted first -0.0613216 0.0918373 -0.668
## hive4 0.0949225 0.0975640 0.973
## IT 0.0588142 0.1421196 0.414
## treatmentWeighted -1.4175699 0.3083761 -4.597
## IT:treatmentWeighted 0.2825651 0.0775090 3.646
##
## Correlation of Fixed Effects:
## (Intr) vstNm_ trt2Wf hive4 IT trtmnW
## vstNm_cntrd -0.055
## trt2Wghtdfr -0.229 0.031
## hive4 -0.038 -0.016 -0.013
## IT -0.993 0.056 0.162 -0.014
## trtmntWghtd -0.312 0.236 0.014 -0.009 0.316
## IT:trtmntWg 0.312 -0.253 -0.013 0.009 -0.318 -0.996ma2 <- update(ma1, .~. - trt2)
BIC(ma1, ma2) # keep ma2## df BIC
## ma1 9 4591.618
## ma2 8 4584.295summary(ma2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log(amp_acc) ~ visitNum_centered + hive + IT + treatment + (1 |
## beeID) + IT:treatment
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 4538.2 4584.3 -2261.1 4522.2 2352
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0293 -0.5416 0.1160 0.6854 2.6651
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.06752 0.2598
## Residual 0.38294 0.6188
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.227873 0.561839 5.745
## visitNum_centered 0.001640 0.000628 2.611
## hive4 0.094133 0.098224 0.958
## IT 0.074144 0.141118 0.525
## treatmentWeighted -1.414972 0.308358 -4.589
## IT:treatmentWeighted 0.281940 0.077505 3.638
##
## Correlation of Fixed Effects:
## (Intr) vstNm_ hive4 IT trtmnW
## vstNm_cntrd -0.049
## hive4 -0.042 -0.016
## IT -0.995 0.052 -0.012
## trtmntWghtd -0.315 0.235 -0.008 0.316
## IT:trtmntWg 0.316 -0.252 0.009 -0.318 -0.996ma3 <- update(ma2, .~. - hive)
anova(ma2, ma3) # keep ma3## Data: newdf3
## Models:
## ma3: log(amp_acc) ~ visitNum_centered + IT + treatment + (1 | beeID) +
## ma3: IT:treatment
## ma2: log(amp_acc) ~ visitNum_centered + hive + IT + treatment + (1 |
## ma2: beeID) + IT:treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma3 7 4537.1 4577.4 -2261.5 4523.1
## ma2 8 4538.2 4584.3 -2261.1 4522.2 0.9083 1 0.3406summary(ma3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## log(amp_acc) ~ visitNum_centered + IT + treatment + (1 | beeID) +
## IT:treatment
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 4537.1 4577.4 -2261.5 4523.1 2353
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0319 -0.5419 0.1149 0.6887 2.6771
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.06929 0.2632
## Residual 0.38295 0.6188
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.250421 0.567438 5.728
## visitNum_centered 0.001648 0.000628 2.625
## IT 0.075722 0.142637 0.531
## treatmentWeighted -1.412882 0.308373 -4.582
## IT:treatmentWeighted 0.281383 0.077509 3.630
##
## Correlation of Fixed Effects:
## (Intr) vstNm_ IT trtmnW
## vstNm_cntrd -0.049
## IT -0.996 0.051
## trtmntWghtd -0.312 0.235 0.312
## IT:trtmntWg 0.313 -0.252 -0.315 -0.996ma3.1 <- update(ma3, .~. - visitNum_centered)
anova(ma3, ma3.1) # keep 3.1## Data: newdf3
## Models:
## ma3.1: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## ma3: log(amp_acc) ~ visitNum_centered + IT + treatment + (1 | beeID) +
## ma3: IT:treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma3.1 6 4541.9 4576.5 -2265.0 4529.9
## ma3 7 4537.1 4577.4 -2261.5 4523.1 6.877 1 0.008731 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(ma3.1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 4541.9 4576.5 -2265.0 4529.9 2354
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0202 -0.5393 0.1093 0.6834 2.6666
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.07006 0.2647
## Residual 0.38402 0.6197
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.32371 0.56952 5.836
## IT 0.05660 0.14315 0.395
## treatmentWeighted -1.60343 0.30014 -5.342
## IT:treatmentWeighted 0.33276 0.07511 4.431
##
## Correlation of Fixed Effects:
## (Intr) IT trtmnW
## IT -0.996
## trtmntWghtd -0.309 0.308
## IT:trtmntWg 0.310 -0.312 -0.996ma4 <- update(ma3.1, .~. - visitNum_centered )
anova(ma3.1, ma4) #keep ma3.1## Data: newdf3
## Models:
## ma3.1: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## ma4: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma3.1 6 4541.9 4576.5 -2265 4529.9
## ma4 6 4541.9 4576.5 -2265 4529.9 0 0 1summary(ma3.1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 4541.9 4576.5 -2265.0 4529.9 2354
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0202 -0.5393 0.1093 0.6834 2.6666
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.07006 0.2647
## Residual 0.38402 0.6197
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.32371 0.56952 5.836
## IT 0.05660 0.14315 0.395
## treatmentWeighted -1.60343 0.30014 -5.342
## IT:treatmentWeighted 0.33276 0.07511 4.431
##
## Correlation of Fixed Effects:
## (Intr) IT trtmnW
## IT -0.996
## trtmntWghtd -0.309 0.308
## IT:trtmntWg 0.310 -0.312 -0.996m3.2 <- update(ma3.1, .~ . + I(scale(visitNum_centered)) + I(scale(visitNum_centered^2)) + (1|beeID))
summary(m3.2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## log(amp_acc) ~ IT + treatment + (1 | beeID) + I(scale(visitNum_centered)) +
## I(scale(visitNum_centered^2)) + IT:treatment
## Data: newdf3
##
## AIC BIC logLik deviance df.resid
## 4539.1 4585.2 -2261.5 4523.1 2352
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0308 -0.5413 0.1154 0.6887 2.6758
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.06929 0.2632
## Residual 0.38294 0.6188
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.26035 0.56770 5.743
## IT 0.07536 0.14273 0.528
## treatmentWeighted -1.41125 0.30921 -4.564
## I(scale(visitNum_centered)) 0.03784 0.01705 2.219
## I(scale(visitNum_centered^2)) 0.00129 0.01797 0.072
## IT:treatmentWeighted 0.28096 0.07774 3.614
##
## Correlation of Fixed Effects:
## (Intr) IT trtmnW I((N_) I((N_^
## IT -0.996
## trtmntWghtd -0.307 0.309
## I(scl(vN_)) -0.057 0.062 0.164
## I(sc(N_^2)) 0.038 -0.035 0.074 -0.511
## IT:trtmntWg 0.308 -0.311 -0.996 -0.177 -0.076anova(ma3.1, m3.2) # keep ma3.1 (according to BIC)## Data: newdf3
## Models:
## ma3.1: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## m3.2: log(amp_acc) ~ IT + treatment + (1 | beeID) + I(scale(visitNum_centered)) +
## m3.2: I(scale(visitNum_centered^2)) + IT:treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma3.1 6 4541.9 4576.5 -2265.0 4529.9
## m3.2 8 4539.1 4585.2 -2261.5 4523.1 6.8821 2 0.03203 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1finmod <- update(ma3.1, .~., REML = TRUE)
summary(finmod)## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## Data: newdf3
##
## REML criterion at convergence: 4545.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0179 -0.5395 0.1071 0.6826 2.6697
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.07463 0.2732
## Residual 0.38434 0.6200
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.32401 0.58488 5.683
## IT 0.05653 0.14700 0.385
## treatmentWeighted -1.60422 0.30031 -5.342
## IT:treatmentWeighted 0.33294 0.07515 4.430
##
## Correlation of Fixed Effects:
## (Intr) IT trtmnW
## IT -0.996
## trtmntWghtd -0.301 0.300
## IT:trtmntWg 0.302 -0.303 -0.996# diagnostics
plot(finmod)
qqnorm(ranef(finmod)$beeID[[1]])
qqline(ranef(finmod)$beeID[[1]])
sjp.lmer(finmod) # plot random effects to find any outliers## Plotting random effects...
sjp.lmer(finmod, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects## `sjp.lmer()` will become deprecated in the future. Please use `plot_model()` instead.## Computing p-values via Wald-statistics approximation (treating t as Wald z).
Bootstrap CI’s for amplitude for figure for paper
# set number of bootstrap replicates for models
nbootSims2 = 1000
# using hive 3, since it's the one with the most data
pframe <- expand.grid(trt2 = levels(droplevels(newdf3$trt2)),
IT = ITmean,
visitNum_centered = -50:50,
hive = factor(3, levels = levels(newdf3$hive)),
beeID = 99999,
treatment = levels(droplevels(newdf3$treatment)))
pframe$amp_acc <- 0
pframe2 <- pframe[(pframe$trt2 == "Sham first" &
pframe$treatment == "Sham" &
pframe$visitNum_centered < 0) |
(pframe$trt2 == "Sham first" &
pframe$treatment == "Weighted" &
pframe$visitNum_centered > 0) |
(pframe$trt2 == "Weighted first" &
pframe$treatment == "Weighted" &
pframe$visitNum_centered < 0) |
(pframe$trt2 == "Weighted first" &
pframe$treatment == "Sham" &
pframe$visitNum_centered > 0)
, ]
# exponentiate to put on original scale
pp <- exp(predict(finmod, newdata = pframe2, re.form=NA, type = 'response')) # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(finmod, FUN=function(x) predict(x, pframe2, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))## [1] "Number of bootstrap samples 1000"bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe2$blo<- exp(bb2_se[1,]) # exponentiate to put on original scale
pframe2$bhi<- exp(bb2_se[2,])
pframe2$predMean <- pp
pframe2## trt2 IT visitNum_centered hive beeID treatment amp_acc
## 1 Sham first 3.965833 -50 3 99999 Sham 0
## 3 Sham first 3.965833 -49 3 99999 Sham 0
## 5 Sham first 3.965833 -48 3 99999 Sham 0
## 7 Sham first 3.965833 -47 3 99999 Sham 0
## 9 Sham first 3.965833 -46 3 99999 Sham 0
## 11 Sham first 3.965833 -45 3 99999 Sham 0
## 13 Sham first 3.965833 -44 3 99999 Sham 0
## 15 Sham first 3.965833 -43 3 99999 Sham 0
## 17 Sham first 3.965833 -42 3 99999 Sham 0
## 19 Sham first 3.965833 -41 3 99999 Sham 0
## 21 Sham first 3.965833 -40 3 99999 Sham 0
## 23 Sham first 3.965833 -39 3 99999 Sham 0
## 25 Sham first 3.965833 -38 3 99999 Sham 0
## 27 Sham first 3.965833 -37 3 99999 Sham 0
## 29 Sham first 3.965833 -36 3 99999 Sham 0
## 31 Sham first 3.965833 -35 3 99999 Sham 0
## 33 Sham first 3.965833 -34 3 99999 Sham 0
## 35 Sham first 3.965833 -33 3 99999 Sham 0
## 37 Sham first 3.965833 -32 3 99999 Sham 0
## 39 Sham first 3.965833 -31 3 99999 Sham 0
## 41 Sham first 3.965833 -30 3 99999 Sham 0
## 43 Sham first 3.965833 -29 3 99999 Sham 0
## 45 Sham first 3.965833 -28 3 99999 Sham 0
## 47 Sham first 3.965833 -27 3 99999 Sham 0
## 49 Sham first 3.965833 -26 3 99999 Sham 0
## 51 Sham first 3.965833 -25 3 99999 Sham 0
## 53 Sham first 3.965833 -24 3 99999 Sham 0
## 55 Sham first 3.965833 -23 3 99999 Sham 0
## 57 Sham first 3.965833 -22 3 99999 Sham 0
## 59 Sham first 3.965833 -21 3 99999 Sham 0
## 61 Sham first 3.965833 -20 3 99999 Sham 0
## 63 Sham first 3.965833 -19 3 99999 Sham 0
## 65 Sham first 3.965833 -18 3 99999 Sham 0
## 67 Sham first 3.965833 -17 3 99999 Sham 0
## 69 Sham first 3.965833 -16 3 99999 Sham 0
## 71 Sham first 3.965833 -15 3 99999 Sham 0
## 73 Sham first 3.965833 -14 3 99999 Sham 0
## 75 Sham first 3.965833 -13 3 99999 Sham 0
## 77 Sham first 3.965833 -12 3 99999 Sham 0
## 79 Sham first 3.965833 -11 3 99999 Sham 0
## 81 Sham first 3.965833 -10 3 99999 Sham 0
## 83 Sham first 3.965833 -9 3 99999 Sham 0
## 85 Sham first 3.965833 -8 3 99999 Sham 0
## 87 Sham first 3.965833 -7 3 99999 Sham 0
## 89 Sham first 3.965833 -6 3 99999 Sham 0
## 91 Sham first 3.965833 -5 3 99999 Sham 0
## 93 Sham first 3.965833 -4 3 99999 Sham 0
## 95 Sham first 3.965833 -3 3 99999 Sham 0
## 97 Sham first 3.965833 -2 3 99999 Sham 0
## 99 Sham first 3.965833 -1 3 99999 Sham 0
## 104 Weighted first 3.965833 1 3 99999 Sham 0
## 106 Weighted first 3.965833 2 3 99999 Sham 0
## 108 Weighted first 3.965833 3 3 99999 Sham 0
## 110 Weighted first 3.965833 4 3 99999 Sham 0
## 112 Weighted first 3.965833 5 3 99999 Sham 0
## 114 Weighted first 3.965833 6 3 99999 Sham 0
## 116 Weighted first 3.965833 7 3 99999 Sham 0
## 118 Weighted first 3.965833 8 3 99999 Sham 0
## 120 Weighted first 3.965833 9 3 99999 Sham 0
## 122 Weighted first 3.965833 10 3 99999 Sham 0
## 124 Weighted first 3.965833 11 3 99999 Sham 0
## 126 Weighted first 3.965833 12 3 99999 Sham 0
## 128 Weighted first 3.965833 13 3 99999 Sham 0
## 130 Weighted first 3.965833 14 3 99999 Sham 0
## 132 Weighted first 3.965833 15 3 99999 Sham 0
## 134 Weighted first 3.965833 16 3 99999 Sham 0
## 136 Weighted first 3.965833 17 3 99999 Sham 0
## 138 Weighted first 3.965833 18 3 99999 Sham 0
## 140 Weighted first 3.965833 19 3 99999 Sham 0
## 142 Weighted first 3.965833 20 3 99999 Sham 0
## 144 Weighted first 3.965833 21 3 99999 Sham 0
## 146 Weighted first 3.965833 22 3 99999 Sham 0
## 148 Weighted first 3.965833 23 3 99999 Sham 0
## 150 Weighted first 3.965833 24 3 99999 Sham 0
## 152 Weighted first 3.965833 25 3 99999 Sham 0
## 154 Weighted first 3.965833 26 3 99999 Sham 0
## 156 Weighted first 3.965833 27 3 99999 Sham 0
## 158 Weighted first 3.965833 28 3 99999 Sham 0
## 160 Weighted first 3.965833 29 3 99999 Sham 0
## 162 Weighted first 3.965833 30 3 99999 Sham 0
## 164 Weighted first 3.965833 31 3 99999 Sham 0
## 166 Weighted first 3.965833 32 3 99999 Sham 0
## 168 Weighted first 3.965833 33 3 99999 Sham 0
## 170 Weighted first 3.965833 34 3 99999 Sham 0
## 172 Weighted first 3.965833 35 3 99999 Sham 0
## 174 Weighted first 3.965833 36 3 99999 Sham 0
## 176 Weighted first 3.965833 37 3 99999 Sham 0
## 178 Weighted first 3.965833 38 3 99999 Sham 0
## 180 Weighted first 3.965833 39 3 99999 Sham 0
## 182 Weighted first 3.965833 40 3 99999 Sham 0
## 184 Weighted first 3.965833 41 3 99999 Sham 0
## 186 Weighted first 3.965833 42 3 99999 Sham 0
## 188 Weighted first 3.965833 43 3 99999 Sham 0
## 190 Weighted first 3.965833 44 3 99999 Sham 0
## 192 Weighted first 3.965833 45 3 99999 Sham 0
## 194 Weighted first 3.965833 46 3 99999 Sham 0
## 196 Weighted first 3.965833 47 3 99999 Sham 0
## 198 Weighted first 3.965833 48 3 99999 Sham 0
## 200 Weighted first 3.965833 49 3 99999 Sham 0
## 202 Weighted first 3.965833 50 3 99999 Sham 0
## 204 Weighted first 3.965833 -50 3 99999 Weighted 0
## 206 Weighted first 3.965833 -49 3 99999 Weighted 0
## 208 Weighted first 3.965833 -48 3 99999 Weighted 0
## 210 Weighted first 3.965833 -47 3 99999 Weighted 0
## 212 Weighted first 3.965833 -46 3 99999 Weighted 0
## 214 Weighted first 3.965833 -45 3 99999 Weighted 0
## 216 Weighted first 3.965833 -44 3 99999 Weighted 0
## 218 Weighted first 3.965833 -43 3 99999 Weighted 0
## 220 Weighted first 3.965833 -42 3 99999 Weighted 0
## 222 Weighted first 3.965833 -41 3 99999 Weighted 0
## 224 Weighted first 3.965833 -40 3 99999 Weighted 0
## 226 Weighted first 3.965833 -39 3 99999 Weighted 0
## 228 Weighted first 3.965833 -38 3 99999 Weighted 0
## 230 Weighted first 3.965833 -37 3 99999 Weighted 0
## 232 Weighted first 3.965833 -36 3 99999 Weighted 0
## 234 Weighted first 3.965833 -35 3 99999 Weighted 0
## 236 Weighted first 3.965833 -34 3 99999 Weighted 0
## 238 Weighted first 3.965833 -33 3 99999 Weighted 0
## 240 Weighted first 3.965833 -32 3 99999 Weighted 0
## 242 Weighted first 3.965833 -31 3 99999 Weighted 0
## 244 Weighted first 3.965833 -30 3 99999 Weighted 0
## 246 Weighted first 3.965833 -29 3 99999 Weighted 0
## 248 Weighted first 3.965833 -28 3 99999 Weighted 0
## 250 Weighted first 3.965833 -27 3 99999 Weighted 0
## 252 Weighted first 3.965833 -26 3 99999 Weighted 0
## 254 Weighted first 3.965833 -25 3 99999 Weighted 0
## 256 Weighted first 3.965833 -24 3 99999 Weighted 0
## 258 Weighted first 3.965833 -23 3 99999 Weighted 0
## 260 Weighted first 3.965833 -22 3 99999 Weighted 0
## 262 Weighted first 3.965833 -21 3 99999 Weighted 0
## 264 Weighted first 3.965833 -20 3 99999 Weighted 0
## 266 Weighted first 3.965833 -19 3 99999 Weighted 0
## 268 Weighted first 3.965833 -18 3 99999 Weighted 0
## 270 Weighted first 3.965833 -17 3 99999 Weighted 0
## 272 Weighted first 3.965833 -16 3 99999 Weighted 0
## 274 Weighted first 3.965833 -15 3 99999 Weighted 0
## 276 Weighted first 3.965833 -14 3 99999 Weighted 0
## 278 Weighted first 3.965833 -13 3 99999 Weighted 0
## 280 Weighted first 3.965833 -12 3 99999 Weighted 0
## 282 Weighted first 3.965833 -11 3 99999 Weighted 0
## 284 Weighted first 3.965833 -10 3 99999 Weighted 0
## 286 Weighted first 3.965833 -9 3 99999 Weighted 0
## 288 Weighted first 3.965833 -8 3 99999 Weighted 0
## 290 Weighted first 3.965833 -7 3 99999 Weighted 0
## 292 Weighted first 3.965833 -6 3 99999 Weighted 0
## 294 Weighted first 3.965833 -5 3 99999 Weighted 0
## 296 Weighted first 3.965833 -4 3 99999 Weighted 0
## 298 Weighted first 3.965833 -3 3 99999 Weighted 0
## 300 Weighted first 3.965833 -2 3 99999 Weighted 0
## 302 Weighted first 3.965833 -1 3 99999 Weighted 0
## 305 Sham first 3.965833 1 3 99999 Weighted 0
## 307 Sham first 3.965833 2 3 99999 Weighted 0
## 309 Sham first 3.965833 3 3 99999 Weighted 0
## 311 Sham first 3.965833 4 3 99999 Weighted 0
## 313 Sham first 3.965833 5 3 99999 Weighted 0
## 315 Sham first 3.965833 6 3 99999 Weighted 0
## 317 Sham first 3.965833 7 3 99999 Weighted 0
## 319 Sham first 3.965833 8 3 99999 Weighted 0
## 321 Sham first 3.965833 9 3 99999 Weighted 0
## 323 Sham first 3.965833 10 3 99999 Weighted 0
## 325 Sham first 3.965833 11 3 99999 Weighted 0
## 327 Sham first 3.965833 12 3 99999 Weighted 0
## 329 Sham first 3.965833 13 3 99999 Weighted 0
## 331 Sham first 3.965833 14 3 99999 Weighted 0
## 333 Sham first 3.965833 15 3 99999 Weighted 0
## 335 Sham first 3.965833 16 3 99999 Weighted 0
## 337 Sham first 3.965833 17 3 99999 Weighted 0
## 339 Sham first 3.965833 18 3 99999 Weighted 0
## 341 Sham first 3.965833 19 3 99999 Weighted 0
## 343 Sham first 3.965833 20 3 99999 Weighted 0
## 345 Sham first 3.965833 21 3 99999 Weighted 0
## 347 Sham first 3.965833 22 3 99999 Weighted 0
## 349 Sham first 3.965833 23 3 99999 Weighted 0
## 351 Sham first 3.965833 24 3 99999 Weighted 0
## 353 Sham first 3.965833 25 3 99999 Weighted 0
## 355 Sham first 3.965833 26 3 99999 Weighted 0
## 357 Sham first 3.965833 27 3 99999 Weighted 0
## 359 Sham first 3.965833 28 3 99999 Weighted 0
## 361 Sham first 3.965833 29 3 99999 Weighted 0
## 363 Sham first 3.965833 30 3 99999 Weighted 0
## 365 Sham first 3.965833 31 3 99999 Weighted 0
## 367 Sham first 3.965833 32 3 99999 Weighted 0
## 369 Sham first 3.965833 33 3 99999 Weighted 0
## 371 Sham first 3.965833 34 3 99999 Weighted 0
## 373 Sham first 3.965833 35 3 99999 Weighted 0
## 375 Sham first 3.965833 36 3 99999 Weighted 0
## 377 Sham first 3.965833 37 3 99999 Weighted 0
## 379 Sham first 3.965833 38 3 99999 Weighted 0
## 381 Sham first 3.965833 39 3 99999 Weighted 0
## 383 Sham first 3.965833 40 3 99999 Weighted 0
## 385 Sham first 3.965833 41 3 99999 Weighted 0
## 387 Sham first 3.965833 42 3 99999 Weighted 0
## 389 Sham first 3.965833 43 3 99999 Weighted 0
## 391 Sham first 3.965833 44 3 99999 Weighted 0
## 393 Sham first 3.965833 45 3 99999 Weighted 0
## 395 Sham first 3.965833 46 3 99999 Weighted 0
## 397 Sham first 3.965833 47 3 99999 Weighted 0
## 399 Sham first 3.965833 48 3 99999 Weighted 0
## 401 Sham first 3.965833 49 3 99999 Weighted 0
## 403 Sham first 3.965833 50 3 99999 Weighted 0
## blo bhi predMean
## 1 27.88884 35.12006 31.19904
## 3 27.97216 35.17425 31.26938
## 5 28.04596 35.23476 31.33989
## 7 28.12669 35.29929 31.41055
## 9 28.20839 35.33684 31.48137
## 11 28.29033 35.39071 31.55236
## 13 28.37251 35.44539 31.62350
## 15 28.45503 35.50759 31.69480
## 17 28.53821 35.57252 31.76626
## 19 28.59370 35.63709 31.83789
## 21 28.68210 35.70151 31.90967
## 23 28.75050 35.76053 31.98162
## 25 28.81907 35.81630 32.05373
## 27 28.88780 35.86353 32.12601
## 29 28.95670 35.90168 32.19844
## 31 29.02589 35.93960 32.27104
## 33 29.11343 35.99425 32.34380
## 35 29.20606 36.06262 32.41673
## 37 29.29943 36.10986 32.48982
## 39 29.39076 36.16208 32.56308
## 41 29.48507 36.21724 32.63650
## 43 29.56848 36.27248 32.71008
## 45 29.64726 36.32781 32.78384
## 47 29.72670 36.38323 32.85776
## 49 29.80635 36.45501 32.93184
## 51 29.87136 36.54930 33.00609
## 53 29.91921 36.64020 33.08051
## 55 29.99823 36.72050 33.15510
## 57 30.09295 36.77759 33.22986
## 59 30.18797 36.83477 33.30478
## 61 30.28329 36.89200 33.37988
## 63 30.35511 36.99573 33.45514
## 65 30.40840 37.11562 33.53057
## 67 30.46235 37.17469 33.60617
## 69 30.52589 37.20508 33.68195
## 71 30.56766 37.29571 33.75789
## 73 30.61995 37.37031 33.83401
## 75 30.67267 37.41339 33.91029
## 77 30.72607 37.47907 33.98675
## 79 30.80508 37.55789 34.06338
## 81 30.91853 37.62946 34.14019
## 83 31.03241 37.69277 34.21716
## 85 31.14670 37.77141 34.29432
## 87 31.23071 37.86450 34.37164
## 89 31.27387 37.95638 34.44914
## 91 31.36409 38.02615 34.52681
## 93 31.43114 38.11712 34.60466
## 95 31.51526 38.19815 34.68269
## 97 31.62292 38.27973 34.76089
## 99 31.71499 38.36068 34.83926
## 104 31.81108 38.52180 34.99655
## 106 31.83995 38.60391 35.07546
## 108 31.90001 38.70058 35.15454
## 110 31.96790 38.78534 35.23381
## 112 32.03549 38.85674 35.31325
## 114 32.10339 38.94356 35.39287
## 116 32.18359 39.04564 35.47267
## 118 32.23909 39.14751 35.55265
## 120 32.29222 39.22410 35.63282
## 122 32.33937 39.30174 35.71316
## 124 32.39874 39.39584 35.79368
## 126 32.47772 39.51021 35.87439
## 128 32.55723 39.60449 35.95527
## 130 32.63758 39.69793 36.03634
## 132 32.72030 39.79161 36.11760
## 134 32.79744 39.86463 36.19903
## 136 32.87484 39.95446 36.28065
## 138 32.93179 40.05022 36.36245
## 140 32.99830 40.13817 36.44444
## 142 33.08852 40.25472 36.52662
## 144 33.18291 40.37152 36.60897
## 146 33.23159 40.48940 36.69152
## 148 33.31052 40.59253 36.77425
## 150 33.37477 40.69400 36.85716
## 152 33.43175 40.79572 36.94027
## 154 33.49066 40.89770 37.02356
## 156 33.56110 40.99993 37.10704
## 158 33.62460 41.10234 37.19070
## 160 33.68271 41.23674 37.27456
## 162 33.73211 41.39638 37.35860
## 164 33.77496 41.52361 37.44284
## 166 33.83217 41.63195 37.52726
## 168 33.88948 41.73962 37.61187
## 170 33.94689 41.84628 37.69668
## 172 34.00439 41.95694 37.78167
## 174 34.06199 42.06651 37.86686
## 176 34.11969 42.17577 37.95224
## 178 34.17749 42.28531 38.03781
## 180 34.23539 42.40321 38.12358
## 182 34.28500 42.52645 38.20954
## 184 34.32780 42.65014 38.29569
## 186 34.37644 42.78285 38.38204
## 188 34.41201 42.89541 38.46858
## 190 34.45351 43.00827 38.55531
## 192 34.52197 43.15890 38.64225
## 194 34.56409 43.31342 38.72938
## 196 34.62097 43.45801 38.81670
## 198 34.66255 43.61892 38.90422
## 200 34.70441 43.77043 38.99194
## 202 34.74647 43.91361 39.07986
## 204 20.55672 26.07949 23.17974
## 206 20.61690 26.11577 23.23201
## 208 20.67727 26.15091 23.28439
## 210 20.73780 26.17571 23.33689
## 212 20.79077 26.20054 23.38951
## 214 20.85701 26.23645 23.44224
## 216 20.92717 26.27437 23.49510
## 218 20.99775 26.31200 23.54808
## 220 21.06808 26.34779 23.60117
## 222 21.13205 26.38790 23.65439
## 224 21.18395 26.42474 23.70772
## 226 21.23607 26.46435 23.76117
## 228 21.28905 26.50667 23.81475
## 230 21.36198 26.54996 23.86845
## 232 21.42087 26.59180 23.92226
## 234 21.47554 26.63316 23.97620
## 236 21.53767 26.67064 24.03026
## 238 21.60035 26.70818 24.08444
## 240 21.66341 26.76016 24.13875
## 242 21.72261 26.80296 24.19317
## 244 21.77216 26.84536 24.24772
## 246 21.83306 26.88782 24.30239
## 248 21.90161 26.93232 24.35719
## 250 21.94748 26.99432 24.41211
## 252 22.00056 27.04006 24.46715
## 254 22.06506 27.07071 24.52232
## 256 22.12992 27.10140 24.57761
## 258 22.19517 27.13248 24.63303
## 260 22.21556 27.16865 24.68857
## 262 22.30524 27.21512 24.74423
## 264 22.39241 27.26625 24.80003
## 266 22.45861 27.31336 24.85594
## 268 22.52501 27.36152 24.91199
## 270 22.57469 27.41286 24.96816
## 272 22.61547 27.46429 25.02445
## 274 22.67089 27.51562 25.08088
## 276 22.73343 27.56726 25.13743
## 278 22.77901 27.61896 25.19411
## 280 22.84112 27.64968 25.25091
## 282 22.88881 27.72290 25.30785
## 284 22.93692 27.77492 25.36491
## 286 22.98559 27.82703 25.42210
## 288 23.03671 27.87924 25.47942
## 290 23.08128 27.94271 25.53687
## 292 23.13520 27.99616 25.59445
## 294 23.18190 28.05627 25.65216
## 296 23.22688 28.12403 25.71000
## 298 23.27534 28.16709 25.76797
## 300 23.32393 28.22074 25.82607
## 302 23.37327 28.26137 25.88430
## 305 23.48898 28.43993 26.00115
## 307 23.53400 28.49247 26.05978
## 309 23.57911 28.54509 26.11854
## 311 23.62463 28.61448 26.17743
## 313 23.69228 28.69139 26.23645
## 315 23.76021 28.76851 26.29561
## 317 23.82834 28.81175 26.35490
## 319 23.89666 28.86456 26.41432
## 321 23.96518 28.93212 26.47388
## 323 24.01281 29.01442 26.53357
## 325 24.07782 29.06770 26.59339
## 327 24.12995 29.14005 26.65336
## 329 24.17227 29.21513 26.71345
## 331 24.22606 29.28164 26.77368
## 333 24.26679 29.35011 26.83405
## 335 24.30759 29.42384 26.89456
## 337 24.35299 29.48817 26.95520
## 339 24.41215 29.56041 27.01597
## 341 24.47146 29.63282 27.07689
## 343 24.51992 29.68743 27.13794
## 345 24.56797 29.75723 27.19913
## 347 24.61660 29.82718 27.26045
## 349 24.66580 29.90754 27.32192
## 351 24.73277 29.96894 27.38352
## 353 24.81203 30.04055 27.44526
## 355 24.85470 30.14112 27.50715
## 357 24.89734 30.21195 27.56917
## 359 24.94006 30.25096 27.63133
## 361 24.96859 30.32426 27.69363
## 363 25.00721 30.39774 27.75607
## 365 25.05565 30.47804 27.81866
## 367 25.10524 30.58094 27.88138
## 369 25.15427 30.66954 27.94424
## 371 25.20482 30.76536 28.00725
## 373 25.24056 30.83883 28.07040
## 375 25.27284 30.89840 28.13369
## 377 25.32440 30.95795 28.19713
## 379 25.37020 31.01827 28.26070
## 381 25.41419 31.09500 28.32442
## 383 25.45712 31.18639 28.38829
## 385 25.50723 31.26283 28.45230
## 387 25.56073 31.31911 28.51645
## 389 25.61107 31.42121 28.58074
## 391 25.65226 31.52449 28.64519
## 393 25.68123 31.62820 28.70977
## 395 25.73948 31.73151 28.77451
## 397 25.78224 31.80616 28.83939
## 399 25.82513 31.87966 28.90441
## 401 25.86815 31.97789 28.96958
## 403 25.91125 32.06483 29.03490### Plot amplitude lmer with trial number
pframe2$trt3 <- mapvalues(pframe2$trt2, from = c("Sham first", "Weighted first"),
to = c("Sham -> Weighted", "Weighted -> Sham"))
g0 <- ggplot(pframe2, aes(x=visitNum_centered, y=predMean))+
geom_line(aes(color = treatment))+
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Sonication number\n (0 is when treatment switched)") +
geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi, fill= treatment), alpha = 0.2) +
facet_wrap(~trt3) +
scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
scale_fill_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
theme(legend.position = c(0.5, 0.89),
legend.background = element_rect(fill=alpha('gray95', 1)),
legend.direction="horizontal")
g0
ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_amp.pdf"), width =5, height = 3.6)
# I like this plot
g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
aes(x = visitNum_centered, y = amp_acc), position = position_jitter(width = 0.5), size = 0.1, pch = 16, alpha = 0.4) +
ylim(c(0, 165))
g1
ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_withRawData.pdf"), width =5, height = 3.6)
# I like this plot
g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
aes(x = visitNum_centered, y = amp_acc), position = position_jitter(width = 0.5), size = 0.1, pch = 16, alpha = 0.4) +
scale_y_log10(limits = c(2, 500)) +
annotation_logticks(sides = "l", color = 'grey', size = 0.2)
g1
ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_withRawData_logScale.pdf"), width =5, height = 3.6)
# make figure that isn't faceted by order of treatment
g03 <- ggplot(pframe2, aes(x=visitNum_centered, y=predMean))+
geom_line(aes(color = treatment))+
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Sonication number\n (0 is when treatment switched)") +
geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi, fill= treatment), alpha = 0.2) +
facet_wrap(~treatment) +
scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
scale_fill_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
theme(legend.position = c(0.75, 0.8),
legend.background = element_rect(fill=alpha('gray95', 1)),
legend.direction="vertical")
g03
ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_trtFacet.pdf"), width =5, height = 3.6)
# I like this plotShow amplitude plot without visit number
look at visit number 1 (right at the time of switching)
# mean and 95% CI (without visit number)
g0 <- ggplot(pframe2[pframe2$visitNum_centered %in% c(1), ], aes(x=treatment, y=predMean))+
geom_point(aes(color = treatment), alpha = 1) +
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment") +
geom_errorbar(aes(x = treatment, ymin = blo, ymax = bhi, color= treatment), alpha = 1, width = 0.1) +
scale_color_viridis(name = "Flower treatment", discrete = TRUE, option = "A", begin =0, end = 0) +
theme(legend.position = "none",
legend.background = element_rect(fill=alpha('gray95', 1)),
legend.direction="horizontal")
g0
ggsave(filename = file.path(figDir, "WeightedSham_Amp_pointWhiskers_BW.pdf"), width =5, height = 3.6)
g0 <- ggplot(pframe2[pframe2$visitNum_centered %in% c(1), ], aes(x=treatment, y=predMean))+
geom_point(aes(color = treatment), alpha = 1) +
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment") +
geom_errorbar(aes(x = treatment, ymin = blo, ymax = bhi, color= treatment), alpha = 1, width = 0.1) +
scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) +
theme(legend.position = "none",
legend.background = element_rect(fill=alpha('gray95', 1)),
legend.direction="horizontal")
g0
ggsave(filename = file.path(figDir, "WeightedSham_Amp_pointWhiskers.pdf"), width =5, height = 3.6)
g22 <- g0 +
annotate(geom="text", x=c(1,2), y=c(0, 0) + 40, label=c("a", "b"),
color="black")
ggsave(filename = file.path(figDir, "WeightedSham_Amp_pointWhiskers_annot.pdf"), width =5, height = 3.6)
## add raw data
g1 <- g0 + geom_point(data = newdf3, aes(x = treatment, y = amp_acc, color = treatment), position = position_jitter(width = 0.1, height = 6), alpha = 0.4, pch = 16, size = 0.2)
g1
ggsave(filename = file.path(figDir, "WeightedSham_amp_pointWhiskers_rawData.pdf"), width =5, height = 3.6)major break – below this line is the first method I used
this is almost identical to the code above
Modeling frequency with lmer
Use BIC to select model (decide what interactions and covariates to use)
Make sure that REML = FALSE when comparing BIC values
Start with the biggest model of interest, and then see what predictors can be removed
# interaction that may be important, based on domain knowledge
m0 = lmer(freq ~ treatment*IT+ hive + (1|beeID), data = new_df, REML = FALSE)
# no interaction
m1 = lmer(freq ~ treatment+IT+ hive + (1|beeID), data = new_df, REML = FALSE)
BIC(m0, m1) # stay with m1 -- no interaction## df BIC
## m0 7 25187.77
## m1 6 25184.69anova(m0, m1) # this p-value is very borderline -- we're being ## Data: new_df
## Models:
## m1: freq ~ treatment + IT + hive + (1 | beeID)
## m0: freq ~ treatment * IT + hive + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 6 25150 25185 -12569 25138
## m0 7 25147 25188 -12567 25133 4.6848 1 0.03043 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# conservative to avoid overfitting by using BIC
m2 = lmer(freq ~ hive + IT+ (1|beeID), data = new_df, REML = FALSE)
BIC(m1, m2) # treatment doesn't make the model better## df BIC
## m1 6 25184.69
## m2 5 25178.53anova(m1, m2) # p-value for paper -- probably don't need to report## Data: new_df
## Models:
## m2: freq ~ hive + IT + (1 | beeID)
## m1: freq ~ treatment + IT + hive + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 5 25150 25178 -12570 25140
## m1 6 25150 25185 -12569 25138 1.6039 1 0.2054# p-value for hive
m3 = lmer(freq ~ treatment + (1|beeID), data = new_df, REML = FALSE)
anova(m1, m3)## Data: new_df
## Models:
## m3: freq ~ treatment + (1 | beeID)
## m1: freq ~ treatment + IT + hive + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 4 25146 25169 -12569 25138
## m1 6 25150 25185 -12569 25138 0.1592 2 0.9235# diagnostics -- use REML = TRUE
m1 <- update(m1, .~., REML =TRUE)
summary(m1) # summary for paper## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ treatment + IT + hive + (1 | beeID)
## Data: new_df
##
## REML criterion at convergence: 25115.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3676 -0.4684 0.1873 0.6825 2.3319
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 954.3 30.89
## Residual 2357.8 48.56
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 355.173 61.937 5.734
## treatmentWeighted 2.614 2.065 1.266
## IT -1.833 15.543 -0.118
## hive4 -4.129 11.398 -0.362
##
## Correlation of Fixed Effects:
## (Intr) trtmnW IT
## trtmntWghtd -0.002
## IT -0.995 -0.016
## hive4 -0.047 0.000 -0.009plot(m1)
qqnorm(ranef(m1)$beeID[[1]])
qqline(ranef(m1)$beeID[[1]])
sjp.lmer(m1, type = "re", sort = TRUE) # plot random effects to find any outliers## Plotting random effects...
sjp.lmer(m1, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects## Computing p-values via Wald-statistics approximation (treating t as Wald z).
Bootstrap CI’s for figure for paper
# set number of bootstrap replicates for models
nbootSims = 1000
table(new_df$hive) # more trials from hive 3##
## 3 4
## 1644 716# using hive 3, since it's the one with the most data
# however, hive doesn't affect model anyway
# calculate an average IT for prediction
ITmean = mean(tapply(new_df$IT, INDEX = new_df$beeID, FUN = function(x) x[1] ))
pframe <- data.frame(treatment = levels(droplevels(new_df$treatment)),
IT = ITmean,
hive = factor(3, levels = levels(new_df$hive)),
beeID = 99999)
pframe$freq <- 0
pp <- predict(m1, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m1, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))## [1] "Number of bootstrap samples 1000"bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe <- pframe[, c('treatment', "blo", "bhi", "predMean")]
pframe## treatment blo bhi predMean
## 1 Sham 335.3492 359.8953 347.9021
## 2 Weighted 337.7346 363.1810 350.5156Make frequency plots for paper
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
g0 <- ggplot(pframe, aes(x=treatment, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Flower treatment") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)
g0
ggsave(plot = g0, filename = file.path(figDir, "HeavyLight_PredsAndCI_freq.pdf"), width =4, height = 3)Modeling amplitude with lmer
# log transform for acceleration helps model fit better
# start with large model (including interaction)
ma0 = lmer(log(amp_acc) ~ treatment * IT + hive + (1|beeID), data = new_df, REML = FALSE)
ma1 = lmer(log(amp_acc) ~ treatment + IT + hive + (1|beeID), data = new_df, REML = FALSE)
BIC(ma0, ma1) # interaction, ma0, is better## df BIC
## ma0 7 4583.334
## ma1 6 4595.150# pval for interaction (in case paper needs it)
anova(ma0, ma1)## Data: new_df
## Models:
## ma1: log(amp_acc) ~ treatment + IT + hive + (1 | beeID)
## ma0: log(amp_acc) ~ treatment * IT + hive + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma1 6 4560.6 4595.1 -2274.3 4548.6
## ma0 7 4543.0 4583.3 -2264.5 4529.0 19.582 1 9.639e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ma2 <- update(ma0, .~. - hive)
BIC(ma0, ma2) # hive not important## df BIC
## ma0 7 4583.334
## ma2 6 4576.547#pval for hive
anova(ma0, ma2)## Data: new_df
## Models:
## ma2: log(amp_acc) ~ treatment + IT + (1 | beeID) + treatment:IT
## ma0: log(amp_acc) ~ treatment * IT + hive + (1 | beeID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma2 6 4541.9 4576.5 -2265.0 4529.9
## ma0 7 4543.0 4583.3 -2264.5 4529.0 0.979 1 0.3224# fit with reml = TRUE for summarizing
ma3 <- update(ma2,.~., REML = TRUE)
summary(ma3) # final model to report## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ treatment + IT + (1 | beeID) + treatment:IT
## Data: new_df
##
## REML criterion at convergence: 4545.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0179 -0.5395 0.1071 0.6826 2.6697
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeID (Intercept) 0.07463 0.2732
## Residual 0.38434 0.6200
## Number of obs: 2360, groups: beeID, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.32401 0.58488 5.683
## treatmentWeighted -1.60422 0.30031 -5.342
## IT 0.05653 0.14700 0.385
## treatmentWeighted:IT 0.33294 0.07515 4.430
##
## Correlation of Fixed Effects:
## (Intr) trtmnW IT
## trtmntWghtd -0.301
## IT -0.996 0.300
## trtmntWg:IT 0.302 -0.996 -0.303# diagnostics
plot(ma3)
qqnorm(ranef(ma3)$beeID[[1]])
qqline(ranef(ma3)$beeID[[1]])
sjp.lmer(ma3) # plot random effects to find any outliers## Plotting random effects...
sjp.lmer(ma3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects## `sjp.lmer()` will become deprecated in the future. Please use `plot_model()` instead.## Computing p-values via Wald-statistics approximation (treating t as Wald z).
Bootstrap CI’s for amplitude for figure for paper
# set number of bootstrap replicates for models
nbootSims2 = 1000
# using hive 3, since it's the one with the most data
pframe <- data.frame(treatment = levels(droplevels(new_df$treatment)), IT = ITmean, hive = factor(3, levels = levels(new_df$hive)), beeID = 99999)
pframe$amp_acc <- 0
# exponentiate to put on original scale
pp <- exp(predict(ma3, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(ma3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))## [1] "Number of bootstrap samples 1000"bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<- exp(bb2_se[1,]) # exponentiate to put on original scale
pframe$bhi<- exp(bb2_se[2,])
pframe$predMean <- pp
pframe <- pframe[, c('treatment', "blo", "bhi", "predMean")]
pframe## treatment blo bhi predMean
## 1 Sham 31.55900 37.97060 34.75029
## 2 Weighted 23.82485 28.84456 26.16290Make amplitude plots for paper
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
# Holding IT = meanIT
print(ITmean)## [1] 3.965833ga0 <- ggplot(pframe, aes(x=treatment, y=predMean))+
geom_point()+
labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
annotate(geom="text", x=c(1,2), y=c(0, 0) + 40, label=c("a", "b"),
color="black")
ga0
ggsave(plot = ga0, filename = file.path(figDir, "HeavyLight_PredsAndCI_amp.pdf"), width =4, height = 3)sessionInfo()## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] viridis_0.4.0 viridisLite_0.2.0 gamm4_0.2-5
## [4] mgcv_1.8-20 nlme_3.1-131 effects_4.0-0
## [7] carData_3.0-0 plyr_1.8.4 multcomp_1.4-8
## [10] TH.data_1.0-8 MASS_7.3-47 survival_2.41-3
## [13] mvtnorm_1.0-6 sjPlot_2.4.0 lme4_1.1-14
## [16] Matrix_1.2-11 reshape2_1.4.2 ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] pbkrtest_0.4-7 RColorBrewer_1.1-2 rprojroot_1.2
## [4] tools_3.4.2 TMB_1.7.11 backports_1.1.1
## [7] R6_2.2.2 sjlabelled_1.0.5 DT_0.2
## [10] lazyeval_0.2.1 colorspace_1.3-2 nnet_7.3-12
## [13] gridExtra_2.3 tidyselect_0.2.3 mnormt_1.5-5
## [16] compiler_3.4.2 sandwich_2.4-0 labeling_0.3
## [19] scales_0.5.0 lmtest_0.9-35 psych_1.7.8
## [22] blme_1.0-4 stringr_1.2.0 digest_0.6.12
## [25] foreign_0.8-69 minqa_1.2.4 rmarkdown_1.8
## [28] stringdist_0.9.4.6 pkgconfig_2.0.1 htmltools_0.3.6
## [31] pwr_1.2-1 htmlwidgets_0.9 rlang_0.1.4
## [34] shiny_1.0.5 bindr_0.1 zoo_1.8-0
## [37] dplyr_0.7.4 magrittr_1.5 modeltools_0.2-21
## [40] bayesplot_1.4.0 Rcpp_0.12.13 munsell_0.4.3
## [43] abind_1.4-5 prediction_0.2.0 stringi_1.1.6
## [46] yaml_2.1.14 merTools_0.3.0 snakecase_0.5.1
## [49] grid_3.4.2 parallel_3.4.2 sjmisc_2.6.2
## [52] forcats_0.2.0 lattice_0.20-35 ggeffects_0.2.2
## [55] haven_1.1.0 splines_3.4.2 sjstats_0.12.0
## [58] knitr_1.17 codetools_0.2-15 stats4_3.4.2
## [61] glue_1.2.0 evaluate_0.10.1 modelr_0.1.1
## [64] nloptr_1.0.4 httpuv_1.3.5 gtable_0.2.0
## [67] purrr_0.2.4 tidyr_0.7.2 assertthat_0.2.0
## [70] mime_0.5 coin_1.2-1 xtable_1.8-2
## [73] broom_0.4.2 survey_3.32-1 coda_0.19-1
## [76] tibble_1.3.4 arm_1.9-3 glmmTMB_0.1.4
## [79] bindrcpp_0.2